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The  purpose  of  this  investigation  was  to  develop  and  verify  a  generalized 
computer  mode]  for  POGO  instability.  The  solution  of  POGO  and  other 
self-excited  instability  problems  is  of  paramount  concern  in  the  design  of  any 
liquid  rocket  system.  Furthermore,  improved  liquid  rocket  systems  will 
almost  certainly  be  a  major  development  item  for  the  support  of  future  Air 
Force  space  systems,  as  well  as  for  the  NASA  space  station. 

A  non-linear  approach  employing  the  method  of  characteristics  in 
conjunction  with  component  models  from  numerous  sources  was  used  in 
generating  a  digital  computer  code  to  determine  both  the  transient  and  steady 
state  responses  of  a  typical  bipropellant  liquid  rocket  to  sinusoidal  thrust 
oscillations.  Since  the  goal  of  this  study  was  to  create  the  computer  code 
and  not  to  perform  an  exhaustive  parametric  analysis,  the  programs  were 
written  in  the  BASIC  language  and  implemented  on  a  personal  computer. 
However,  this  code  could  easily  be  transported  to  a  main  frame  system  (such 
as  a  CYBER)  and  a  thorough  design  study  would  be  most  enlightening. 
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Notation 


A: 

Element  Cross  Sectional  Area 

a: 

Sonic  Velocity 

A': 

Dorsch  Pump  Parameter 

®inj* 

Dorsch  Injector  Constant 

BO: 

Pump  Coefficient 

Bl: 

Pump  Coefficient 

B2: 

Pump  Coefficient 

B3: 

Pump  Coefficient 

B*: 

Dorsch  Pump  Parameter 

C*: 

Effective  Exhaust  Velocity 

c  *■ 

S  • 

Nominal  Exhaust  Velocity 

Cslopo 

*:Local  Slope  of  C 

CM: 

Backward  Characteristic  Value 

cp: 

Forward  Characteristic  Value 

CT; 

Thrust  Coefficient 

C*: 

Dorsch  Pump  Parameter 

D: 

Effective  Pipe  Diameter 

Valve  Closure  Parameter 

F: 

Thrust 

■v 

Nominal  Thrust 

f 

Friction  Factor 

g 

Local  Acceleration 

i 

max 

Number  of  Single  Pipe  Sections 

K: 

Fluid  Bulk  Modulus 

L*: 

Chamber  Characteristic  Length 

M: 

Mass 

M: 

Mass  Flow  Rate 

MR: 

Mixture  Ratio 

NS: 

Number  of  Sections  in  Element 

P: 

Pressure 

P  - 

Combustor  Pressure 

PPo‘- 

Nominal  Pump  Inlet  Pressure 

R: 

Gas  Constant  (exhaust  gas) 

Tc: 

Combustor  Temperature 

*cuf 

Valve  Closure  Time 

V: 

Velocity 

Voi: 

Cavitation  ‘Bubble’  Volume 

Vo,0: 

Initial  'Bubble'  Volume 

V 

Nominal  Pump  Flow  Velocity 

X: 

Displacement 

§:  Area  Ratio 

r 

Pipe  Material  Young's  Modulus 

AHjnj:  Dorsch  Injector  Pressure  Drop 

AP: 

Pump  Pressure  Rise 

APq:  Nominal  Pump  Pressure  Rise 

At: 

Time  Increment 

Ax:  Section  Separation 

0: 

Grid-Mesh  Ratio 

0 Gas  Residence  Time 

9 

P: 

Density 

T:  Orifice  Flow  Constant 

Tc: 

Chamber  Dead  Time 

fa):  Circular  Frequency 

C:  Grid  Parameter 

Un: 

Natural  Frequency 

Subscripts 


c: 

Combustor 

o: 

Downstream 

Dischg:  Discharge  Line 

F: 

Fuel  System 

Fm£ 

Feed  Line 

inj: 

Injector 

Ok: 

Oxidizer  System 

p: 

Pump 

T: 

Tank 

u: 

Upstream 

uii: 

Ullage 

Pressure  and  Velocity  Indices 

Pressure  and  velocity  terms  for  the  systems  modeled  in  this  study  are 
designated  with  indices  ana  qualifying  subscripts.  The  indicies,  I  thru  7 
refer  to  the  Figs.  4.  IB,  4.3B,  and4.6B  and  the  single  pipe,  monopropellant, 
and  bipropellant  systems,  respectively.  Qualifying  identifiers  are  found  in 
the  subscripts  list. 
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In  this  thesis,  a  generalized  computer  code/for  the  simulation  of  POGO 

instability  in  liquid  rockets  was  generated  anti  verified.  The  term  POGO 

/ 

refers  to  the  out  of  phase  motion  of  the  ends  of  a  liquid  rocket,  which  can 

build  to  a  dangerous  magnitude  due  to  propulsion  feed  back.  The  elimination 

i 

of  any  POGO  instability  will  be  of  great  importance  in  the  development  of 
future  liquid  rocket  systems  such  as  QTV^rdntH&W* .  The  model  detailed  in 
this  investigation  is  based  on  a  method  of  characteristics  solution  of  the 
simplified  Navier-Stokes  equations.  The  resulting  non-linear  differential 
equations  were  solved  using  first  order  finite  difference  methods.  This 
solution  was  applied  to  the  fluid  lines  of  several  liquid  rocket  systems. 
Boundary  conditions,  based  on  component  models  used  generally  throughout 
the  literature  were  developed  to  link  these  fluid  lines. 

The  computer  routines  were  verified  by  comparison  with  published 
results  from  several  sources.  The  results  of  numerous  runs  agreed  quite 
well  with  the  published  data,  even  in  very  non-linear  systems  with  large 
disturbance  amplitudes. 

As  an  aid  to  future  investigators  the  routines  developed  in  this  thesis 
were  applied  to  a  typical  bipropellant  liquid  rocket  system.  -"'Both  the 
ransient  and  the  steady  state  responses  of  this  system  to  sinusoidal  thrust 
ons  were  obtained.  A  relatively  limited  parametric  study  was 
perf«>. .  and  indicated  that  this  particular  system  was  very  stable  and 

showed  no  sign  of  POGO  instability.  Two  factors  were  most  Important  in  the 
system's  stability:  the  pump  operating  points  and  the  chamber  dead  time;  the 
stability  of  the  system  being  greatly  enhanced  by  the  the  stable  range  in 
which  the  pumps  operated  and  the  selection  of  a  relatively  short  combustor 
dead  time. 

A  user's  guide  was  compiled  to  detail  the  software  developed  in  this 
investigation.  Its  purpose  is  to  facilitate  the  application  of  these  routines  to 
other  systems  by  future  investigators.  _ 
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Solving  for  VD  in  terms  of  quadratic  coeficients, 

For  VD: 

B  =  x2g(pau-aD)  [3.8] 

C  =  t2(Cm-Cp);  p  = 

Au 

The  quadratic  coefficients  A,  B,  and  C  are  used  with  [3.6]  in  determining 
V.  The  other  unknowns  may  now  be  computed.  This  boundary  condition  is 
similar  to  the  section  area  change  except  that  a  specified  loss  is  implied  by 
Eq.  [3.7]. 

Relative  Component  notion.  Fig.  3.3  depicts  relative  motion  between 
two  system  elements.  If  the  areas  and  orientations  of  the  two  elements  are 
dissimilar,  there  will  in  general  be  a  finite  rate  of  fluid  storage  at  their 
interface  due  to  the  relative  motion.  Equating  the  rate  of  storage  to  the  rate 
of  inflow  minus  the  rate  of  outflow  at  the  interface1*'^ 


Figure  3.3 

Relative  component  motion. 

vu  Au~  vd  ad  =  -VsCAySinau-  ADsinaD)  [3. 9] 
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Using  the  simple  continuity  relation  for  an  incompressible  fluid14: 


Vu=advd  (3-31 

Assuming  no  energy  dissipation  (i.e.  assuming  no  losses)  at  the 
interface14137  and  equating  the  total  pressure  at  either  side  of  the  interface: 

pu*^evu2  =  p0*'/2pvD2  [3.4] 

There  are  now  four  equations  in  four  unknowns.  Solving  in  terms  of  the 
downstream  velocity: 

.  AVd2+BVd+C  =  0  [3.5] 

where:  p  =  ^2  A  =  (1-p2) 

Au 

B  =  2<vP«u>  C  =  |(C„-Cp) 

The  unsubscripted  terms,  A,  B,  and  C  are  quadratic  coefficients.  From 
the  sense  of  Fig.  3.Z,  A  will  be  a  positive  quantity.  Thus  to  obtain  a  positive 
flow  velocity,  the  positive  square  root  should  be  applied  in  the  quadratic 
formula :  _ 

v-  -B»v/B2-4Ac'  [3.61 

2  A 

Note  that  other  boundary  conditions  will  be  described  in  terms  of 
quadratic  coefficients  (the  identifiers  A,  B,  and  C  will  be  used  in  each  case). 
Unless  otherwise  stated,  use  of  the  positive  form  of  the  quadratic  equation 
([3.6])  will  be  assumed.  It  is  now  possible  to  solve  for  the  other  quantities. 
This  completely  specifies  conditions  at  both  end  points.  Note  that  the  only 
assumption  made  was  to  neglect  energy  dissipation  at  the  interface.  For  the 
typical  liquid  rocket,  this  is  a  reasonable  assumption  due  to  the  relatively 
small  velocities  in  the  feed  system  as  well  as  the  streamlined  geometries 
involved1*"3. 


Known  Pressure  Upstream  fir  Downstream.  If  the  pressure  at  the 
upstream  or  downstream  end  of  a  conduit  is  a  known  function  of  time  F(t), 
this  relation  can  be  combined  with  Eqs.  [2.16]  or  [2.17]  to  fix  the  endpoint 
conditions14.  For  the  upstream  and  downstream  cases 


A 


Known  Pressure  Upstream 

f[j=F(t);  =  fy ;  $=  CM+paDVD 

Known  Pressure  Downstream 

^  =  F(t);  ly=  ly  =  Cp-payVy 


13.1] 


Note  that  the  subscripts  (J  and  D  refer  to  points  fnfinttesstmally  upstream 
and  downstream  of  the  end  in  question,  not  to  tne  end  itself.  Combining  the 
equations  to  find  the  velocity  at  either  end  of  the  pipe: 


F(t)-  CM 

V-T5T3  ,3’n 


V 


CD-F(t) 

P»u 


[3.2] 


Eqs.  [3.1]  and  [3.21  give  the  velocities  at  either  end  of  the  pipe.  Since 
the  pressure  has  already  been  specified,  the  boundary  conditions  are 
completely  fixed.  Both  the  upstream  and  downstream  elements  In  this 
system  will  affect  each  other.  Thus,  the  dynamic  effects  of  each  element  in 
a  complex  system  will  propagate  both  upstream  and  downstream. 

Section  Area  Change.  Fig.  3.2  depicts  the  interface  between  lines  of 
constant  area.  Each  line  is  modeled  as  in  the  previous  section  and  it  is 
necessary  to  develop  the  relations  which  connect  them.  Applying  a 
compatibility  equation  at  each  pipe  end 


Figure  3.2 


Section  area  change. 


III.  Elements  and  Boundary  Conditions 


in  this  section,  the  major  components  of  a  typical  liquid  rocket  will  be 
examined  with  the  aim  of  applying  the  knowledge  of  their  dynamic  behavior 
(taken  from  various  sources  in  the  literature)  to  the  problem  of  mating  them 
to  their  connecting  fluid  lines.  First  however,  it  will  be  necessary  to 
describe  the  connection  between  these  external  elements  and  the  finite 
difference  equations  which  model  the  pipe  flow  portion  of  the  problem.  That 
is,  a  general  form  for  the  various  types  of  boundary  conditions  in  the  system 
must  be  determined. 


Boundary  Conditions 

At  either  end  of  a  single  pipe  only  one  of  the  compatibility  equations 
(Eqs.  [2.16]  or  [2.17])  is  available.  This  situation  is  depicted  in  Fig. 
3.1 1^.36.  At  the  upstream  end  of  the  pipe,  only  the  backward  difference 
equation  (Eq.  [2.17])  is  available;  the  situation  is  reversed  at  the 
downstream  end.  To  determine  the  pressure  and  velocity  at  the  ends,  it  will 
be  necessary  to  develop  auxiliary  equations  (boundary  conditions)  based  on 
end  conditions. 


To  maintain  convergence,  theCourant  condition14  must  be  satisfied,  namely 


0  -  — —  C  £  ~2— 

°  V+a'  ^  V+a 


[2.27] 


This  simply  Implies  that  in  Fig.  2. 1 ,  the  points  R  and  S  fait  between  A  and 
B.  Finally,  combining  Eqs.  [2. 16]  and  [2. 17]  and  solving  for  the  pressure  at 
point  P14: 

_  CP+CM 

=  12.28] 


Either  of  Eqs.  [2.16]  or  [2. 17]  can  be  used  to  compute  the  velocity.  This 
completely  establishes  the  state  at  all  internal  points  in  the  pipe.  Note  that  a 
linear  interpolation  has  been  used.  In  order  to  maintain  the  accuracy  of  the 

non-linear  results,  the  values  of  £  and  0  should  approach  the  Courant  limits, 
implying  only  a  small  interpolation.  The  pipe  flow  problem  is  completely 
solved  at  this  point,  however  the  problem  of  establishing  the  boundary 
conditions  at  the  end  points  (where  either  Cp  or  Cm  are  unknown)  remains. 


importance,  at  Boundary  Contitloiw 

The  equations  developed  in  this  section  allow  the  computation  of  the 
pressure  and  velocity  at  essentially  arbitrary  points  inside  a  conduit. 
Unfortunately,  a  liquid  rocket  is  more  than  a  simple  collection  of  pipes. 
Indeed,  the  boundary  conditions  make  up  the  most  active  components  of  the 
system.  In  the  next  section,  boundary  conditions  which  link  the  pipe 
elements  of  a  rocket  propulsion  system  will  be  developed. 


highly  deformable  tubes  or  compressible  fluids,  the  variation  in  At  would 
greatly  complicate  a  complex  system  with  multiple  elements  since  boundary 
conditions  must  be  passed  from  one  element  to  the  next.  For  this  reason 
constant  time  steps  will  be  used. 

Specified  Time  Intervals.  Rearranging  Eqs.  [2. 12J  and  [2.  H]14: 


Noting  that 


Similarly, 


^=Cp-paRVp 

12.161 

Pp=CM-pasV0 

[2.17] 

Cp=P^eV4'4Atsina'  fA2DVRl  ] 

R 

[2.18J 

r  Q  fAt|Vcl  ■ 

Cm-  *s  +  P®svs[l+v“  ^tsina  2d 

s 

|  [2.19] 

,  1  It  is  apparent  that 

XC-*R  VVR 
*c-x*'Vc-VA 

[2.20] 

VC-<R(VC-VA) 

R'  1  *aR<vc-vA> 

[2.21] 

8=5*  12221  4=SSa  =  0. 

12.23] 

vc-Cs(vc-vy 

s"  1  +  0s(Vc-Vb) 

[2.24] 

Pr=Pc-(VrBr.«r)(Pc-Pa) 

[2.25] 

Ps=Pc  +  (Vs0s  +  Cs)  (PC-PA) 

[2.26] 

9 
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The  pressure  and  velocity  at  points  A,  C,  and  B  which  occur  at  time  t  are 
known  (either  from  the  last  time-step  iteration  or  from  steady  state 
results).  Conditions  at  points  R  and  S  occur  at  time  t  but  must  be  calculated 

from  known  relationships.  The  state  at  point  P  occurs  at  time  t+At.  Our 
approach  is  to  use  Eqs.  [2.8H2.10]  to  compute  the  conditions  at  P. 
Multiplying  by  the  differential  time  dt,  and  integrating  (to  first  order): 


Figure  2.1 

Twng-displdoffwni  grid. 


*p'fR+  P °r(Vp- VR)  -  p eRg since  (tp-tR) 
JvRIVRl , 

♦poR-fsr-WiO 


12.12] 


Xp-*R=  (VR-aR)(tp-tR) 


12.131 


Pp-Ps-pas(Vp-Vs)  +  pasgs1n<x  (tp-ts) 

fvJVJ, 


(2.14] 


Xp-xs=(Vs-es)(tp-ts) 


[2.15] 


There  are  two  possible  approaches  in  obtaining  a  numerical  solution  to 
these  equations:  the  use  of  a  grid  of  characteristics  or  the  use  of  specified 
time  Intervals.  With  the  grid  of  characteristics,  the  time-step  is  chosen  so 
that  the  points  R  and  S  fall  on  points  A  and  B  respectively.  While  this  can 
lead  to  improved  accuracy  (relative  to  the  use  of  specified  time  intervals)  in 


[2.1]  and  [2.2].  If  the  appropriate  selections  are  made,  simplification  is 
possible.  In  particular,  since  P  and  V  are  functions  of  x  and  t,  regarding  x 
as  a  function  of  t,  then 


—  -P  —  +P  i^-v~+V+ 

rlt  ~  **  rt*  +  '  (Jt  Vx  **  V* 


dt 


dt  -  'x  dt 

This  can  be  accomplished  if 


Now,  Eq.  [2.3]  becomes  the  ordinary  differential  equation 


[2.4! 


[2.51 


Solving  Eq.  [2.5] 


,  dP  dV  fvlvl 

A  dt  dt  2D 


\  -  _] _ .  d*  -  v  -f  a 

pa  '  dt  "V 

A  “  pa  '  dt  “  V 


12.6] 


[2.7] 


Substituting  Eq.  [2.7]  Into  [2.0]  produces  the  following  total  differential 
equations: 

dP  dV 
dt  p  dt 


dP 

dt 


dV 

'dt 


-  pagstna  +  pa 

O 

II 

>a 

>CM 

[2.8] 

d*  _  v 
dt  V  +  a 

[2.9] 

+  pagslnoc  -  pa 

o 

II 

>CN 

[2.10] 

42U  V-a 
dt  V  8 

12.11] 

Finite  Difference  Equations 


Traditionally,  a  finite  difference  approach  is  used  to  solve  these 
decidedly  non-linear  equations.  Fig.  2.1  is  a  time-displacement  grid 

depicting  the  fluid  conditions  at  various  points  in  a  pipe  at  times  t  and  t+At. 


II.  Method  of 


The  method  of  characteristics  transforms  partial  differential  equations 
for  which  no  general  solution  exists  (such  as  those  describing  liquid  flow  in  a 
conduit)  into  particular  total  differential  equations.  The  resulting  non-linear 
equations  may  then  be  integrated  using  numerical  methods  employing  finite 
difference  equations. 

In  this  section,  the  equations  describing  continuity  and  momentum  for 
incompressible  fluids  will  be  transformed  into  total  differential  equations. 
Forward  and  backward  finite  difference  equations  will  then  be  developed  for 
implementation  with  constant  time-step  computations. 


Emaiiops  sL  Motion 

The  hydraulic  equations  embodying  the  principles  of  conservation  of 
momentum  and  continuity  in  a  one  dimensional  pipe  are  respectively14 

5uwx+Vt-gsinoc  ♦  =0  [2.1] 

f»+f>V  +  ga2Vx  =  0  [2.2] 

These  general  equations  can  be  combined  with  an  unknown  multiplier  A, 
to  yield 

*[i(v+ jp)  *p.]  *  r  3] 

-  g  since  ♦  =  0 

Any  choice  of  two  real,  distinct  values  of  A  will  again  yield  two 
independent  equations  in  terms  of  P  and  V  that  are  still  equivalent  to  Eqs. 


► 

Bioropellant  Liquid  Rocket.  The  methods  developed  in  this  investigation- 
(  will  be  applied  to  a  complete,  closed  loop,  liquid  rocket  burning  liquid 

hydrogen  and  oxygen  propellants.  The  system  analyzed  is  intended  to  be  a 
typical  one  and  will  serve  as  a  base  line  for  future  investigators. 

Risulti.  The  data  resulting  from  the  single  pipe  and  monopropellant 
>  computer  models  will  be  compared  with  published  data  to  establish  the 

accuracy  of  the  routines  generated  in  this  thesis.  The  transient  and  steady 
state  response  of  the  bipropellant  rocket  to  various  inputs  will  be  presented. 
An  exhaustive  parametric  analysis  of  these  responses  is  beyond  the  scope  of 
*  this  study. 

Conclusions.  The  overall  perfomance  of  the  models  used  in  this 
investigation  will  be  examined  in  the  light  of  the  objectives  of  this  thesis. 
Future  uses  for  the  routines  developed  in  this  thesis  will  be  recommended. 
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This  thesis  may  best  be  broken  up  into  five  sections.  The  intent  is  to 
provide  both  a  step-by-step  derivation  of  the  system  model  and  a  reference 
for  future  implementations  of  the  programs  developed  here. 

Method  of  Characteristics  Solution.  To  provide  a  framework  for  a 
finite  difference  solution,  the  necessary  equations  to  model  fluid  flow  in  feed 
lines  will  be  derived  using  the  method  of  characteristics14.  These  equations 
form  the  very  heart  of  this  analysis.  With  a  clear  understanding  of  the  pipe 
flow  component  of  this  problem,  the  other  parts  become  simple  boundary 
conditions. 

Boundary  Conditions.  In  order  to  complete  the  feed  line  solutions 
obtained  using  the  method  of  characteristics,  boundary  conditions  must  be 
specified.  These  conditions  take  on  a  definite  physical  significance  in 
describing  the  valves,  pumps,  accumulators,  and  other  components  which 
make  up  a  liquid  rocket.  In  this  study,  we  will  employ  ‘typical*  models, 
i.e.,  element  models  which  have  been  extensively  used  in  the  literature. 
Specific  knowledge  of  a  given  component  is  of  great  importance  in  improving 
the  overall  accuracy  of  any  predictive  tool.  Future  investigators  will 
undoubtedly  wish  to  employ  more  complex  boundary  conditions  to  improve  the 
overall  model  accuracy. 

Systems  Analysis.  In  this  section,  the  results  derived  above  will  be 
used  to  model  several  systems.  In  each  case,  a  system  model  will  bef 
assembled  and  Implemented  in  the  BASIC  computer  language. 

Single-Pine  System.  A  single-pipe  system  will  be  designed  and 
several  different  boundary  conditions  will  be  applied'in  order  to  gain  physical 
insight  into  their  effects.  In  the  next  section,  results  will  be  compared  with 
those  of  Wylie  and  Streeter14,  who  analyzed  an  identical  system. 

Monooropellant  System.  A  complex  open  loop  system  first  analyzed 
by  Oorsch,  Wood  and  Lightner1  will  be  modeled  with  relations  developed  in 
this  paper.  This  model  will  be  used  subsequently  for  verification  by 
comparison  with  published  results. 


propellant  feed  components  including  lines,  valves,  and  orifices  with  the  aim 
of  applying  these  techniques  to  high  frequency  combustion  instability.  He 
showed  that  disturbances  led  to  large  (i.e.  unlimited)  amplitude  waves  in 
closed  pipe  systems  and  to  limited  amplitudes  In  pipe  flow  systems  such  as 
liquid  rocket  motors. 

In  1965,  Fashbaugh  and  Streeter?  used  the  method  of  characteristics  with 
various  boundary  conditions  to  form  one  of  the  first  complete  models  of  a 
propellant  feed  system  in  the  analysis  of  the  observed  POGO  instability  of  the 
Titan  II  missile.  The  authors  used  finite  difference  methods  to  simulate  the 
transients  involved  on  a  digital  computer.  Fashbaugh  and  Streeter  concluded 
with  the  suggestion  that  the  fluid  system  results  could  be  coupled  with  engine 
thrust  relations  and  structural  dynamic  equations  to  complete  the  closed- 
loop  system.  The  major  thrust  of  this  paper  is  to  close  the  loop  and  obtain 
the  overall  transient  and  steady  state  response  of  a  liquid  rocket  to  arbitrary 
thrust  oscillations. 

More  recently,  Dorsch  et  alJ  of  the  NASA  Lewis  Research  Center  have 
employed  the  wave  plan,  a  pulse  synthesis  method,  fundamentally  similar  to 
the  method  characteristics,  to  examine  the  response  of  a  simple 
monopropellant  feed  system.  Boundary  conditions  for  important  feed  system 
components  were  developed  and  the  authors  also  published  a  great  deal  of 
experimental  results  for  a  monopropellant  system. 


QlUgctiyre 

The  purpose  of  this  investigation  is  three  fold:  to  develop  a  closed  loop 
model  of  POGO  instability  from  the  well  established  method  of 
characteristics,  to  verify  the  model  by  comparison  with  results  published  by 
other  analysts,  and  to  provide  a  reference  and  computer  routines  for  future 
investigators  of  similar  systems.  These  objectives  will  be  achieved  only  if 
the  routines  generated  here  are  both  accurate  enough  to  be  a  useful  design 
tool,  and  flexible  enough  to  handle  a  varied  group  of  systems. 


This  can  result  in  a  large  closed  loop  gain  and  the  build-up  of  dangerously 
large  vibrations.  Most  OTV  and  OMV  designs  differ  radically  from  current 
vehicles.  Since  OTV*s  and  OMVs  are  designed  to  function  solely  outside  of  the 
atmosphere,  they  have  a  much  lower  profile,  being  relatively  short  and 
squat.  The  magnitude  of  P060  instability  in  these  vehicles  will  be  of  great 
interest. 


Previous  Work 

The  P060  phenomenon  and  its  associated  complications  have  been  the 
subject  of  considerable  study.  Two  basic  types  of  analysis  have  been 
performed  to  date:  linearized  network  analyses  and  time-step  simulations. 
Linear  network  analyses  employ  perturbation  techniques  on  linearized 
impedances  at  mean  operating  conditions,  while  time-step  methods  employ  a 
non-linear  method  of  characteristics  solution  solved  using  a  finite  difference 
approach. 

Linearized  Models.  Network  methods  of  analysis  have  been  employed 
by  Ryan2,  Rubin3**0,  Holster  and  Astleford3,  and  Zielke6  (superscripts 
designate  references  given  at  the  end  of  the  paper).  These  investigators 
used  linear  mathematical  models  to  establish  transfer  functions  for  stability 
analysis.  Network  methods  are  Ideal  for  the  analysis  of  complex,  but 
primarily  linear,  systems.  After  the  component  transfer  functions  are 
determined,  a  conventional  analysis  (generally  using  the  Nyquist  stability 
criterion)  can  be  performed.  The  simplifications  Involved  are  usually  valid 
for  the  purposes  of  making  closed  loop  stability  predictions  where  the  growth 
or  decay  of  small  sinusoidal  perturbations  is  considered.  Because  of  the 
extremely  nonlinear  nature  of  liquid  rocket  systems,  linearized  models 
cannot  be  used  to  find  the  wave  shapes  of  large  amplitude  pressure  and  flow 
disturbances,  nor  can  they  find  the  response  associated  with  events  such  as 
valve  closures*12. 

Time-Step  Methods.  In  an  early  work,  Woods4  used  a  method  of 
characteristics  formulation  to  model  flow  fluctuations  in  a  number  of 
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I.  Introduction 

Th8  development  of  high  performance  interorbitai  transfer  vehicles  will 
enable  a  major  expansion  of  space  utilization  and  dramatically  reduce  the 
overall  cost  of  space  operations.  To  accomplish  objectives  such  as  these, 
OTYs  (Orbital  Transfer  Vehicles)  and  OMYs  (Orbital  Maneuvering  Vehicles) 
are  being  designed  to  provide  a  new  level  of  performance  in  chemical 
rockets.  These  vehicles  will  support  future  Air  Force  systems  and  NASA's 
space  station  well  into  the  21st  century. 

Propulsion  systems  utilizing  liquid  oxygen  and  hydrogen  offer  significant 
inherent  advantages  over  other  technologies.  They  can  easily  achieve  higher 
performance  (principally  higher  specific  impulse)  and  versatility  than  solid 
propellant  systems.  Furthermore,  they  represent  more  mature  technologies 
than  exotic  systems  such  as  ion,  arcjet  and  plasmajet  concepts.  There  are, 
of  course,  many  technical  difficulties  involved  in  achieving  gains  In 
performance.  As  always,  a  thorough  understanding  of  these  difficulties  is 
prerequisite  to  their  solution. 

A  significant  difficulty  in  the  development  of  any  liquid  rocket  system  is 
the  elimination  of  self-excited  instability.  One  of  the  pricipal  concerns  is  the 
clear  understanding  and  reduction  of  so  called  'POGO*  instabilities.  The  name 
POGO  is  derived  from  the  similarity  observed  between  the  out  of  phase  motion 
of  the  ends  of  a  rocket  vehicle  and  the  motion  of  a  Pogo  stick.  This  vibration 
arises  from  the  interactions  of  the  propellant  feed  system,  the  system  thrust 
function,  and  the  rocket  component  structures.  These  oscillations  have  been 
observed  in  a  number  of  important  launch  vehicles  including  the  Titan  II, 
Thor/Agena,  and  Saturn  SIIB.  During  a  POGO  event,  thrust  variation,  typical 
in  the  normal  operation  of  virtually  all  liquid  rockets,  leads  to  structural 
and  propellant  oscillations  which  further  Intensify  the  original  disturbance. 

1 


I 


The  structural  velocity,  Vs  is  assumed  to  be  directed  vertically,  while 
the  fluid  velocities,  Yy  and  VD  are  directed  longitudinally  with  their 

respective  fluid  lines  (at  orientation  angles  CCy  and  «D).  Thus,  Yg,  Vy,  and 

YD  are  scalars.  Using  Eqs.  12.16]  and  12. 171,  and  noting  that  the  pressure  is 
essentially  constant  across  the  Interface: 

PU=  CP~PauVU;  *D  =  CM+PaDVD 

PU  =  PD 

Solving  these  equations  simultaneously: 


V  PVD_^S(sinctu"  PsinaD) 
where:  r  =  — 


[3.10] 


The  other  unknowns  may  now  be  readily  determined.  This  completes  the 
description  of  the  fluid  dynamics  of  relative  motion. 

Accumulator  in  Lino.  An  accumulator  (as  in  Fig.  3.4)  adds  a 
localized  compliance  to  the  interface  between  two  elements.  Since  an 
accumulator  is  an  energy  storage  device  (the  total  energy  and  fluid  storage 
being  a  function  of  pressure)  the  local  pressure  will  be  described  by  a 
differential  equation.  The  compliance  of  an  accumulator  is  defined  by,:8 


Figure  3.4 

Accumulator  in  lino . 


[3.11] 


where:  a' =  compliance 
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As  in  the  relative  motion  case,  the  fluid  storage  rate  in  the  accumulator 
is  equal  to  the  inflow  rate  minus  the  outflow  rate,:1 


V  A  -  V  A  - 
VU  D  “ 


dVoi 

dt 


[3.121 


Applying  the  calculus  and  substituting 

dVoi_d^  dP  .  dP 

dt  “  dP  dt  dt 

Therefore 

dPp  VuVVpA, 
dt  "  a1 


[3.13] 


This  gives  the  time  history  of  the  pressure  at  the  accumulator  as  a 
function  of  the  flow  rates.  The  pressure  may  be  found  using  any  of  several 
numerical  integration  schemes  (a  simple  first  order  integration  is  used  in 
this  analysis).  Once  the  pressure  is  determined,  the  boundary  condition  is 
exactly  the  same  as  that  for  a  known  pressure. 

Instantaneous  Pressure  Change.  The  action  of  a  pump  may  be 

approximated  by  a  lumped  pressure  rise  of  magnitude  AP  if  the  pump's 
dimensions  are  small  compared  to  the  other  system  elements.  The  resulting 
boundary  condition  is  very  similar  to  the  known  pressure  case: 


ly=  Cp- pflyVy;  =  CM+paDVD 


PD=  I*  +  AP 


Auvu  =  ad^d 


Solving  simultaneously 


Cn-  +  AP 
VD=  P  ---tl 


P(0V°D> 


[3.14] 


Again,  the  other  quantities  can  readily  be  calculated.  Note  that  no  loss 
term  has  been  used  at  the  interface.  Any  actual  losses  should  be  taken  into 

account  when  assigning  a  functional  form  to  AP. 
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System  Components 


A  typical  liquid  rocket  system  consists  of  several  different  components 
which  give  rise  to  the  boundary  conditions  described  above.  The  models 
found  in  this  section  are  used  generally  throughout  the  literature. 

Tanlc  Pressurization.  The  oxidizer  and  fuel  tanks  will  be  pressurized 
to  a  specified  pressure  known  as  the  ullage  pressure.  This  is  done  to  ensure 
that  a  minimum  pressure  (net  positive  suction  head)  is  maintained  in  the  feed 
system  in  order  to  avoid  excessive  pump  inlet  cavitation.  For  purposes  of 
this  analysis,  the  ullage  pressure  will  be  assumed  constant.  Thus  the 
boundary  condition  consists  of  a  simple  known  pressure. 

Valves.  A  valve  can  accurately  be  modeled  as  an  orifice  with  a  variable 

flow  constant  X.  The  boundary  condition  associated  with  the  valve  is 
precisely  the  same  as  that  developed  for  the  orifice. 

Pump  Inlet  Cavitation.  This  phenomenon  is  very  poorly  understood. 
Various  authors  have  used  differing  approaches,  most  simply  adjusting  the 
sound  speed  in  the  pump  inlet  so  that  the  numerical  analysis  matched 
experimentally  observed  resonance  conditions^**?.  Others  have  modeled  a 
hypothetical  accumulator  at  the  pump  inlet  to  account  for  the  locally  high 
compliance  of  the  cavitation  *bubble*,»,°  (this  should  not  be  confused  with 
pump  cavitation  and  the  surge  phenomena  as  a  small  volume  of  vapor  is 
normally  present  in  any  cryogenic  pump  inlet3**).  The  latter  approach  is 
used  in  this  analysis  simply  because  it  seems  less  artificial.  One  must 
recognize  however,  that  this  Is  still  an  approximation  and  the  only  meaningful 
criterion  to  be  used  in  judging  the  success  of  this  model  is  comparison  with 
experimental  results.  The  relative  amount  of  inlet  compliance  must  be 
adjusted  so  that  the  results  are  In  agreement  with  experimental  data. 

The  compliance  associated  with  the  formation  of  vapor  bubbles  (due  to 
their  isothermal  expansion  and  contraction)  in  the  pump  inlet  is  given  by 
Oorsch,  et  al. 1:8  as 


This  relation  need  only  be  substituted  into  Eq.  [3.13]  to  obtain  the 


fc 


boundary  condition 

d*U  _  PU  %  AU~  [3  t  g] 

«  "  PoVol, 

Thus,  the  accumulator  is  positioned  at  the  end  of  the  upstream  pipe. 
Again,  this  highly  non-linear  equation  is  solved  using  a  first  order  numerical 

integration  technique.  Specifcally,  Py  is  integrated  to  first  order  while  the 

right  hand  side  of  [3.16]  is  evaluated  at  mean  conditions  for  the  time  step. 
These  mean  conditions  are  determined  by  iteration  (using  a  modified  form  of 
Newton's  method). 

Inlector  Pome  Compliance.  The  injector  is  housed  in  a  dome  with  a 
compliance  due  to  the  cavitation  of  the  cryogenic  propellant.  For  this 
reason,  Dorsch,  et  al.  modeled  this  compliance  in  exactly  the  same  way  as 
the  pump  inlet  compliance111 1  ( i . e.  with  an  isothermal  vapor  bubble 
expansion  model).  The  same  analysis  will  be  used  in  this  investigation. 

Pump  Pressure  Wise.  The  dynamic  pressure  output  of  a  turbo  pump  is 
generally  similar  to  Figs.  3.3  and  3.6.  These  depict  the  variation  of  pump 
output  pressure  with  output  velocity  (or  flow  for  incompressible  fluids)  and 
suction  pressure,  respectively.  The  normal  operating  point  is  near  the  peak 
in  Fig.  3.5  and  well  away  from  the  cavitation  point  in  Fig.  3.6.  An  excellent 
curve  fit  for  such  pump  characteristics  is  given  by  Fashbaugh  and 
Streeter711014  as 
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Figure  3.5 

Pressure  rise- flow  characteristic. 


Figure  3.6 

Pressure  rise-inlet  pressure  characteristic. 


ap  =  ap0-b0(vd-vd)2-(b,- 


<pu 


Pu.> 


B,-B,(P,  -P,r 


[3.17] 


■U  'U0 

The  coefficients  above  are  determined  from  the  characteristic  plots 
(Figs.  3.5  and  3.6).  The  boundary  condition  associated  with  the  pump  has 
been  derived  previously. 

Combuator.  The  thermochemistry  of  the  combustor  is  relatively  well 
understood.  Several  authors  have  devised  a  composite  time-lag  theory;  a 
complete  description  of  which  is  beyond  the  scope  of  the  present 
work1  18»20s  Dorsch,  et  alJ:**  used  a  two-part  time-lag, 

associating  a  constant  dead  time  with  the  burning  of  the  propellant  and  a 
residence  time  (a  function  of  chamber  geometry,  exhaust  velocity,  and 
temperature)  with  propellant  injection.  The  chamber  pressure  rate  of 
change  is  given  by 


dt  BnVol 


Bf 


i*r* 

[3-'91 

In  applying  these  equations  to  a  bipropellant  system  the  total  injector 
mass  flow  rate  will  be  the  sum  of  the  mass  flow  rates  for  the  fuel  and  the 
oxidizer.  This  differential  equation  establishes  another  boundary  condttlon  to 
be  numerically  integrated.  For  purposes  of  this  analysis  the  effective 


exhaust  velocity,  C*  is  considered  a  linear  function  of  the  instantaneous 
mixture  ratio,  MR*: 

C*  =  C**C*slop#(MR-MR0)  [3.20] 


In  reality,  for  a  bipropellant  liquid  rocket  the  effective  exhaust  velocity 
is  a  very  complex  function  of  mixture  ratio  as  well  as  chamber  pressure  and 
temperature.  This  simplification  Is  justified,  however,  If  the  range  of 
instantaneous  mixture  ratios  is  sufficiently  small  as  will  be  the  case  in  the 
bipropellant  system  analyzed  in  this  investigation. 


Implementation 

The  boundary  conditions  associated  with  the  various  elements  described 
in  this  section  will  be  combined  with  the  fluid  flow  model  of  Section  II  to 
dynamically  model  a  liquid  rocket.  Various  simple  systems  will  first  be 
implemented  to  gain  insight  into  the  behavior  of  each  component  before  they 
are  combined  in  more  complex  arrangements. 
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In  this  section,  the  finite  difference  equations  and  boundary  conditions 
detailed  previously  will  be  used  to  model  three  complete  systems  of 
increasing  complexity.  A  very  simple,  single  pipe  system  will  be  modeled  to 
verify  the  accuracy  of  a  digital  computer  program  written  using  the  equations 
given  in  Sections  it  and  IK  by  comparison  with  published  results.  Various 
boundary  conditions  will  then  be  applied  to  this  simple  model  to  gain  an 
understanding  of  the  effects  of  several  parameters. 

A  more  complex  monopropellant  rocket  system  first  modeled  by  Dorsch, 
et  a).1  will  next  be  examined.  The  POGO  analysis  rouines  developed  in  this 
investigation  will  be  applied  to  this  system  and  results  will  be  generated  for 
comparison  with  those  published  by  Dorsch,  et  al. 

Finally,  to  demonstrate  a  more  general  application  of  the  software 
developed  in  this  investigation  to  systems  of  arbitrary  complexity,  a 
hypothetical  bipropellant  liquid  rocket  system  will  be  designed  and  modeled 
with  these  POGO  analysis  routines.  A  very  limited  parametric  analysis  will 
be  performed  to  determine  those  factors  important  in  POGO  instability  (see 
Section  V). 


Single  Pine  Models 

Fig.  1.1A  depicts  a  simple  pipe  flow  problem  with  a  constant  pressure 
reservoir  upstream  and  a  valve  downstream.  Table  4. 1  lists  values  for  the 
parameters  involved. 


Table  4.1 

Single  Pipe  Parameters 

8=1,200  m/s 

L=600  m 

D=0.5  m 

5=1.4715  MPa 

Tq=0.00205 

f=0.018 

Tcut=2.1  s 

F>=0  MPa 

Ec=!.5 

The  values  in  Table  4.1  were  taken  from  an  identical  system  modeled  by 
Wylie  and  Streeter1*38.  Initially,  the  system  is  assumed  to  have  reached 
steady  state  conditions  (with  the  valve  open;  its  flow  constant  being  taken  as 

Tq).  These  steady  state  conditions  are  easily  calculated  using  the  simplified 

Navier-Stokes  equations  [2.1]  through  [2.2]  since  the  partial  temporal 
derivatives  are  zero  here  by  definition  (see  results  In  Section  V).  At  time 
equal  to  zero,  the  valve  is  shut  according  to  the  following  exponential  law: 


This  system  (in  a  vertical  rather  than  horizontal  orientation  with  respect 
to  the  local  acceleration  vector)  and  two  other  sytems  adding  a  time-varying 
reservoir  pressure  and  variable  acceleration  magnitude  (simulating  the 
effects  of  rocket  thrust  variation)  are  also  examined  in  response  to  a 
sinusoidally  varying  reservoir  pressure.  For  these  responses,  the  orifice 
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flow  constant,  T  is  taken  as  the  constant,  TQ  (implying  that  the  valve 
remains  open). 

Boundary  Conditions,  in  order  to  model  the  system  using  the  finite 
difference  equations  of  Section  II  boundary  conditions  from  Section  III  must 
be  applied.  Fig.  4.  IB  schematically  depicts  the  single  pipe  and  its  two 
interfaces  (with  the  reservoir  and  valve): 


< 

)  1 

Pipe 

Reservoir 

Valve  *, 

Discharge 

Conditions  * 

Conditions  * 

Conditions 

(Index  0,  U) 

(Index  1 ,  U) 

(Index  1 ,  D) 

Figure  4. IB 

Single  pipe  Interfaces. 


Note  that  by  the  index  notation  of  the  figure,  the  reservoir  pressure  and 
velocity  of  the  system  are  denoted  P0>u  and  V0<u,  respectively  while  the  valve 

pressure  is  P,jj  and  the  valve  flow  velocity  is  The  destinction  between 

conditions  just  upstream  and  downstream  of  an  interface  is  necessary  since 
some  of  the  boundary  conditions  in  Section  III  involve  spatial  discontinuities 
In  pressure  and  flow  velocity. 

Boundary  Condition  $£  Interface  fl.  This  is  the  known  upstream 
pressure  condition  described  by  equations  [3.11  and  [3.2].  In  terms  of  Fig. 
4.1  A  indices 


P0,D  =  *0,U 


V0,  D  =  V0,  U  = 


=  Po 

Pq-Cm 

pa 


[4.2] 


Note  that  is  always  available  since  conditions  inside  the  pipe  are 
known  (either  from  the  previous  time  step  or  from  the  steady  state 
solution). 
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Boundary  Condition  at  Interface  1 .  This  condition  is  a  combination  of 
the  known  downstream  pressure  condition  of  [3.1]  and  [3.2]  and  the  orifice 
condition  described  by  [3.7]  (with  flow  constant  modeled  by  [4. 1]).  For  the 
case  at  hand 


For  V, 


u>- 


A  =  1 


B  =  T2pt»p 

c  =  -x2cP  p  =  ^i-  =  i 

A  pip* 

V2 


[4.3] 


,  D=  0 


Note,  V,  u  is  given  in  terms  of  quadratic  coefficients  for  solution  by 

[3.6]. 

Other  Single  Pine  Models.  To  examine  some  of  the  effects  of  boundary 
conditions  more  typical  of  a  liquid  rocket,  another  single  pipe  system  was 
modeled.  The  system  is  based  on  that  taken  from  Wylie  and  Streeter  (Fig. 
4. 1A,B)  but  differs  as  noted  below. 

Single  Pipe  Compliant  System.  For  this  system  the  valve  transient 
response  is  deleted  and  the  reservoir  pressure  is  given  the  following 

sinusoidal  form  (oscillating  at  frequency  u) 


P0=  1.4715  MPa  +100  KPa  -  smut 


[4.4] 


An  accumulator  Is  placed  in  line  with  the  valve  to  model  the  effects  of 
valve  inlet  cavitation.  The  isothermal  expansion  model  (see  [3.15])  was 
used  for  the  compliance.  This  changes  the  boundary  condition  at  interface  I 
to  a  combination  of  the  known  downstream  pressure  ([3.1]  and  [3.2]),  orifice 

([3.7]),  and  the  accumulator  ([3.13])  conditions  (P,  y  becomes  the 

accumulator  pressure  which  must  be  integrated  numerically).  In  terms  of 
the  pressures  and  velocities  of  Fig.  4.  !A: 


[4.5] 


Note  that  the  average  values  which  are  used  in  the  first  order  integration 
of  the  accumulator  pressure  are  determined  by  iteration  (see  Appendix  A  for 

a  source  code  listing  of  this  routine).  P^,  the  initial  bubble  pressure,  was 

taken  as  the  steady  state  pressure  at  the  valve,  143  kPa.  Three  different 
initial  accumulator  volumes  were  examined:  t,  0.1,  and  0.01  m^.  The 
purpose  of  these  changes  is  to  produce  the  kind  of  response  to  pressure 
variations  seen  at  the  injectors  of  a  liquid  rocket  (the  injectors  are  housed  in 
a  dome  with  a  cavitation  compliance). 

Program  Flow.  Any  finite  element  routine  will  involve  a  great  many 
repetitive  calculations,  necessitating  an  organized  approach.  The  flow 
diagram  for  the  problem  at  hand  is  depicted  in  Fig.  4.2: 


Figure  4.2 

Singlt  pip*  computer  flow  diagram. 
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The  program  consists  of  four  main  parts:  initialization,  initial  boundary 
conditions,  interior  sections,  and  terminal  boundary  conditions.  The 
initialization  includes  reading  in  constants  and  steady  state  results  from  a 
simple  static  analysis.  The  pressure  and  velocity  at  the  first  point  in  the 
pipe  are  calculated  using  the  constant  pressure  boundary  condition.  Mote 
that  the  old  pressure  and  velocity  must  be  available  for  the  calculation  of 
conditions  at  the  next  element;  they  are  therefore  held  in  memory.  The  state 
at  each  interior  section  is  calculated  next,  updating  the  old  pressures  and 
velocities  at  each  point.  Finally,  the  state  at  the  last  pipe  element  is 
calculated  using  the  orifice  boundary  condition.  See  the  Appendix  for  a 
source  code  listing. 


Dorsch,  et  al. '  analyzed  a  liquid  rocket  system  using  the  so-called  'wave 
plan'  developed  at  the  NASA  Lewis  Research  Center.  This  pulse  synthesis 
approach  is  closely  related  to  the  method  of  characteristics  used  in  this 
study  and  a  comparison  of  results  from  the  model  of  this  monopropellant 
system  developed  in  this  section  with  those  of  Oorsch,  et  al.  will  verify  the 
accuracy  of  the  routines  detailed  in  Appendix  A.  The  system  is  depicted  in 
Fig.  4.3A. 

The  System.  The  monopropellant  system  consists  of  the  following 
components:  fuel  tank,  feed  line,  turbo  pump,  throttle  valve,  discharge  line, 
injector,  and  combustor.  A  cavitation  bubble  is  hypothesized  to  exist  in  the 
pump  inlet  and  has  a  compliance  modeled  with  the  isothermal  expansion 
results  of  Section  III.  The  system  is  assumed  to  be  at  equilibrium  initially, 
then  it  is  excited  by  the  forced  sinusoidal  motion  of  the  propellant  tank  and 
pump.  As  there  is  no  propulsion  loop  feed  back  and  the  dynamic  nature  of  the 
pump  and  tank  motions  is  ignored,  the  system  may  be  considered  to  be  open. 

Parameter  Values.  Parameter  values  describing  the  monopropellant 
system  are  displayed  in  Table  4.2  .  These  values  were  taken  from  Dorsch,  et 


alJ*,7f  and  represent  a  typical  monopropellant  liquid  rocket.  Note  that  the 
sonic  velocity  is  different  in  each  pipe  due  to  varying  degrees  of  pipe 
elasticity.  The  acoustic  velocity  within  a  pipe  is  given  by*:^63 


[4.7] 


Other  parameters  in  Table  4.2  are  defined  as 

/  C*  L*  \ 

Combustion  Parameter  =  f - )  Au 

VVolcg0g/  AHp 

Pump  Pressure  Rise,  AHp  *  A+BQ  +  CQ^ 

Injection  Velocity, 

Boundary  Conditions.  The  component  interfaces  for  the 
monopropellant  liquid  rocket  which  must  be  modeled  with  boundary  conditions 
are  depicted  in  Fig.  4.3B: 


0  1  2  3  4  5 


Tank 

Feed 

Line 

Feed 

Line 

Discharge 

Line 

Discharge 

Line 

Feed  Pump  Discharge  Injector  & 
_ Valve _ Valve  Combustor 

Figure  4.3B 

Monopropallant  system  component  interfaces. 


Interface  £  (Ullage  Pressure-Propellant  Tank).  This  is  a  known 
upstream  pressure  condition  (the  tank  ullage  pressure  is  assumed  constant) 
modeled  exactly  as  interface  0  of  the  single  pipe  system. 
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Figure  4.10 

Bipropellant  liquid  rocket  system  flow  diagram. 


Figure  4.8 

Pump  pressure  rise  response  to  pump  flow  variation. 


Figure  4.9 

Pump  pressure  rise  response  to  inlet  pressure  variation. 


_ Table  4.4 _ 

_ Bipropellant  Structural  Parameters _ 

Structural  Properties 

Natural  Frequency  10  Hz  Damping  Factor  0.1 

(Undamped ) 

Mass  Spring  Constant  Damping  Constant 
(Kg)  (MN/m)  (KN-sec/m) 

Structure  3000  - 

Tank 

Fuel  3236  12.774 

Ox.  19414  76.643 

Pump 

Fuel  100  0.39478 

Ox.  100  0.39478 


Table  4.5 

Bipropellant  System  Calculated  Parameters 

Char.  Length,  L* 

7  m 

Chamber  Vol 

.  ,  Voir  0 

.  1  m3 

Chamber  Temp. ,  Tc  3500  K 

Gas  Constant 

,  Ra  5.8867  J/Kg  •  K 

Nom.  Accel. ,  g  19. 

Maee  Flow  Rate 

32  m/sec2 

Fuel 

16.18 

Kg/sec 

Oxidizer 

97.07 

Kg/sec 

Injector  Constant 

29.664-103 

(mW'2 

30.787-103 

<m3/Kg>,/2 

Injector  Diameter 

45 

mm 

27 

mm 

Nominal  Pump 

Pressure  Rise,  APq 

58.153 

MPa 

58.006 

MPa 

Nominal  Pump 

Flow  Velocity,  Vp 

12.896 

m/sec 

15.006 

m/sec 

Nominal  Pump 

Inlet  Pressure,  1* 

100.24 

KPa 

326.93 

KPa 

Pipe  Wall  Thickness 

Tank 

2.0448 

mm 

1 . 9245 

mm 

Feed  Line 

0.64437 

mm 

0.30505 

mm 

Oischanrge  Line 

1.9718 

mm 

0.50658 

mm 

Pump  Coefficients 

BO 

50000 

Kg/m2 

50000 

Kg/m2 

B1 

0 

Pa 

0 

Pa 

B2 

0.5 

Pa 

0.5 

Pa 

B3 

5-10" 

6 

5-10" 

6 

40.663 

243.96 

1 .2566 
1.2566 
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Table  4.3 


Bipropellant  Fluid  System  Parameters 


Thrust 


500  KN 


Nominal  Mixture  6:1 

(CT*. /Fuel  Mass  Rates) 


Properties 

Exhaust  Velocity,  C*  4415  m/sec 
(Effective)  at  6M 

4562  m/sec 
at  8:1 


Thrust  Coef.,  Cj  0.999  Chamber  Pressure,  35  MPa 

Residence  Time,  8g  1.5  msec  Dead  Time,  4  msec 


Material  Properties 


Young's  Modulus,  V 
(Pipe  Material ) 


Density,  p 
(Fluid) 

Bulk  Modulus,  K 


10  GPa 

Fuel  (L  H2) 
71  Kg/m3 

02.2  MPa 


Oxidizer  (LOx ) 
1  140  Kg/m3 

1.641  GPa 


Injector  Pres.  Drop 
Ullage  Pressure 


Fluid  system  Properties 

40V.  Injection  Velocity 


150  m/sec 


150  KPa 


Tank 


Length  0iamet< 
(m)  (cm) 

6.448  300 

2.409  300 


Sound  Speed  Friction  Factor  Alpha 


Feed  Line 
Fuel 
Ox. 

Discharge  Line 
Fuel 


( m/sec  ) 


0.03 

0.03 

0.03 

0.03 

0.002 

0.002 


(Deg ) 


39 


described  by  13.20].  The  system  thrust  must  also  be  calculated  and  is  given 
by20:357 


F  = 


cTc#(rt5  jD  F+rt5DjQx) 


[4.13] 


Tank  and  Pump  Motions.  The  motions  of  the  tanks  and  pump  depicted  in 
Fig.  4.7  are  described  by  the  following  total  differential  equations  (these 
follow  directly  from  Newton's  third  law): 


"SXS  =  XT,F  kT,F  +  XT,F  CT,F  +  XT,0x  kT,0x+XT,0xC 
+  XP,F  kP,F  +  XP,F  CP,F  +  XP,0xkP,0x+XP,0xC 

Mt  f(Xs+Xt  f)=-Xt  f  kLF“XT,F  CT,F 
M T, Ox^ XS+  X =  “X T,0xkT , Ox  XT, 0xC T, Ox 

Mp  pCXg*1  Xp  p)  — — Xp  F  kp  p— Xp  p  Cp  p 

MP, oi*S+  *P, oi=_XP, O^P,  Ox  XP, OxCP, Ox 


T,Ox 
P,  Ox 


[4.14] 


It  is  convenient  to  Integrate  these  equations  numerically  at  each  time 
step  (although  several  sub  time  steps  are  used  for  improved  accuracy)  to 
obtain  the  component  velocities  (for  use  in  the  relative  component  motion 
boundary  conditions)  and  the  structure  acceleration,  gv  (for  use  in  the  finite 
difference  equations). 

Program  Flow.  Fig.  4.10  is  the  computer  flow  diagram  for  the 
bipropellant  rocket  POGO  analysis  routine.  It  bears  considerable  similarity 
to  the  monopropellant  case,  except  for  the  thrust  determination  and 
structural  feed  back  routines.  A  user's  guide  to  this  program  is  available  in 
Appendix  A. 
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4.9  depict  the  pump  pressure  rise,  AP,  response  to  pump  flow  velocity. 


V3  0,  and  pump  inlet  pressure,  P3  y,  variations  respectively.  These  curves 

arise  from  the  choice  of  pump  coefficients  for  use  with  [3. 17].  The  rationale 
behind  this  choice  of  coefficients  is  to  provide  a  stable  range  of  pump 
operation  for  this  analysis  (see  Appendix  B). 

Boundary  Conditions.  Fig.  4.6B  defines  the  component  interfaces  for 
the  purpose  of  establishing  boundary  conditions  for  the  bipropellant  POGO 
analysis  routine. 
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BipropelUnt  system  component  interfaces. 


A  comparison  with  Fig.  4.3B  shows  that  the  fluid  elements  of  the 
bipropellant  system  are  exactly  the  same  as  those  of  the  monopropellant 
system  (for  which  boundary  conditions  have  already  been  developed),  except 
that  the  bipropellant  system  contains  both  fuel  and  oxidizer  subsystems.  In 
fact,  the  boundary  conditions  for  the  bipropellant  system  can  be  obtained  by 
simply  adding  the  subscripts  F  and  Ox  to  each  term  employing  indices  from 

Fig.  4.3B  i.e.  P,  y  becomes  P^  y  F  and  P1  y  .  The  combustor  conditions 

are  the  same  except  that  the  total  mass  flow  rate  (fuel  plus  oxidizer)  must 
be  substituted  for  the  monopropellant  mass  flow  rate,  and  the  effective 
exhaust  velocity  is  a  function  of  the  instantantaneous  mixture  ratio  as 
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Structure-component  interactions . 


in  visualizing  the  complex  interactions  between  model  components.  Note  that 
the  thrust  feedback  closes  the  propulsion-structural  vibrations  loop.  The 
only  new  problems  are  the  thrust  calculation  and  dynamic  response  of  the 
floating  (relative  to  the  vehicle  structure)  tanks  and  pumps.  The  system 
thrust  is  related  to  the  total  injection  mass  flow  rate  as  well  as  to  the 
chamber  specific  impulse  (itself  a  function  of  the  instantaneous  mixture 
ratio).  The  thrust  variation  feeds  back  through  the  structure  to  drive  the 
propellant  tank  and  pump  oscillations. 

Model  Components.  Fig.  4.6A  is  a  simplified  view  of  the  bipropellant 
feed  system.  It  is  quite  similar  to  two  monopropetlant  models  in  parallel, 
having  both  fuel  and  oxidizer  components.  Fig.  4.7  depicts  the  bipropellant 
structural  system.  This  is  a  five  degree  of  freedom  model  allowing  pump  and 
tank  motions  relative  to  the  structure.  The  'structure*  (denoted  by  the 
subscript  *s‘)  includes  all  elements  except  the  pumps  and  tanks  (i.e.  all 
other  fluid  lines,  vehicle  structure,  and  any  payload). 

Choice  at  Perimeters.  Before  a  numerical  analysis  may  be  made, 
values  must  be  chosen  for  all  system  parameters.  Tables  4.3  and  4.4  give 
values  chosen  for  this  bipropellant  system.  These  values  represent  a 
'typical*  system  and  do  not  correspond  to  any  physical  implementation. 
Table  4.5  gives  values  for  other  parameters  calculated  from  Tables  4.3  and 
4.4  using  simple  steady  state  relationships 1 the  details  of  this  static 
analysis  as  well  as  the  reasons  for  all  parameter  values  chosen  are 
available  in  Appendix  B.  The  fuel  (liquid  hydrogen)  and  oxidizer  (liquid 
oxygen)  components  of  the  bipropellant  system  are  quite  similar,  except  that 
the  oxidizer  is  more  than  ten  times  as  dense  as  the  fuel^O  and  the  sonic 
velocities  in  the  oxidizer  lines  are  much  less  than  those  for  the  fuel 
component  due  to  the  relative  compliance  (i.e.  the  wall  elasticity  of  the  fluid 
lines  compared  to  the  bulk  modulus  of  the  fluid)  of  the  oxidizer  pipes.  These 
contrasting  values  were  selected  primarily  to  demonstrate  the  effects  of 
sonic  velocity  on  the  performance  of  each  (fuel  and  oxidizer)  fluid  system. 
The  fuel  and  oxidizer  pump  parameters  will  prove  to  be  of  principal 
importance  in  determining  the  systems  over  all  POGO  stability.  Figs.  4.8  and 
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Figure  4.4 

Monopropellant  sysbotn  computer  flow  diagram. 


Bloropfillant  Liquid  Bacfcal 

In  order  to  Illustrate  the  application  of  the  routines  and  principles 
detailed  in  this  thesis,  a  closed  loop  model  for  a  typical  liquid  rocket  was 
developed.  Each  fluid  system  (fuel  and  oxidizer)  of  the  bipropellant  model  is 
similar  to  the  monopropellant  system  of  Dorsch,  et  al.  with  several 
important  differences: 

1.  The  system  thrust  is  calculated  from  chamber  conditions  and  is  fed 
back  into  the  system  resulting  in  a  time-varying  acceleration7. 

2.  The  propellant  tanks  and  pumps  are  modeled  dynamically  as  separate 
mass-spring-dashpot  systems3*4'7.  Their  relative  motions  are  driven  by 
variations  in  the  vehicle  acceleration. 

1.  The  system  is  driven  by  sinusoidal  thrust  variations  imposed  upon 
(simply  added  to)  Item  1  above7. 

Schematic.  The  system  schematic  diagram,  Fig.  4.5,  is  given  as  an  aid 
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Interface  5.  (Discharge  Line-Combustion  Chamfer).  The  combustor 
pressure  is  solved  for  separately  (at  the  beginning  of  each  time  step)  and 
can  be  considered  known,  thus  this  condition  involves  the  orifice  equation, 
[3.7],  a  compatibility  equation,  [2.16],  the  section  area  change  equation, 
[3.3]  (the  discharge  line  area  is  quite  different  than  the  injector  area),  and 
the  accumulator  equation,  [3.13],  to  model  injector  dome  compliance. 
Combining  these  equations  and  integrating  the  injector  inlet  pressure 
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[4.12] 


The  other  unknowns  are  readily  determined.  Note  that  the  average 
combustor  pressure,  Pc  ^  Is  simply  the  mean  of  the  current  and  most 

recent  (the  result  of  the  last  time  step)  combustor  pressures. 

Combustor  Pressure.  The  combustor  pressure  is  numerically 
integrated  using  [3. 18].  This  is  relatively  straight  forward  since  the  injector 
mass  flow  rate  can  be  determined  from  old  values  of  the  injector  velocity, 
V5  D  (see  [3.18]).  These  old  values  are  available  from  a  first-in-first-out 

(FIFO)  queue  which  is  maintained  by  the  routine.  This  queue  holds  a  complete 
set  of  injection  velocities  for  a  time  period  equal  to  the  chamber  dead  time, 

V 

Program  Flow.  Fig.  4.4  Is  a  greatly  simplified  flow  diagram  for  the 
routine  used  to  model  the  monopropellant  rocket  system.  The  fuel  tank,  feed 
line,  and  the  two  discharge  lines  (separated  by  the  throttle  valve),  are 
modeled  as  simple,  one-dimensional  pipes.  The  flow  through  each  of  these 
lines  is  treated  exactly  as  in  the  single  pipe  example,  except  that  the  initial 
boundary  condition  of  one  line  is  the  terminal  boundary  condition  of  the  next. 
The  results  from  this  program  will  be  compared  with  those  of  Dorsch,  et  al. 1 
in  Section  V. 
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compliance  and  relative  motion.  This  involves  the  compatibility  relations, 
[2. 16]  and  [2.17],  the  accumulator  equation,  [3. 13]  (to  model  the  pump  inlet 
cavitation)  modified  by  the  relative  component  motion  equation,  [3.9],  and 
the  instantaneous  pressure  change  equation,  [3.14].  in  terms  of  Fig.  4.3B 
indices,  these  are 
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These  constitute  a  set  of  four  equations  in  four  unknowns.  Solving 
and  integrating  P3  y  to  first  order 
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The  pump  output  pressure,  P3  Q  is  determined  from  [3.14].  The 

average  values  must  be  found  through  an  iterative  process 

Interface  4  (Discharge  Line-Discharge  Line).  This  boundary  condition 
is  exactly  the  same  as  for  Interface  2  with  the  subscript  4  substituted  for 
subscript  2. 
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Interface  1  (Propellant  Tank-Feed  Line).  This  boundary  condition 
involves  the  section  area  change  equations,  [3.2]  through  [3.4]  except  that 
the  simple  incompressible  continuity  relation,  [3.3],  is  replaced  by  the 
relative  component  motion  continuity  equation  ,  [3.9].  Combining  these 


equations  and  solving  for  V,  D  in  terms  of  quadratic  coefficients  (for 


subtitutfon  into  [3.6]) 
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After  V{  Q  has  been  determined,  the  other  unknowns  at  the  boundary 

are  readily  calculated  using  [3.9],  [3.4],  and  a  compatibility  relation  (either 
[2.16]  or  [2.17]). 

Interface  2.  (Feed  Line-Feed  Line).  This  boundary  condition  Is 
modeled  as  the  simple  orifice  condition  of  [3.7]  and  [3.8]  with  an  orifice 

constant  equal  to  Tv^  (the  feed  valve  flow  constant).  In  terms  of  Fig.  4.3B 


For  V2  D: 
A  =  1 


p 

af*«j 


B  =  Ty,  P  P(°2,  U  ~  Pa2,  ^  =  -"^V1  Cp) 


[4.9] 


a  Again,  V2  0  Is  found  by  substitution  into  [3.6]  and  the  other  unknowns 

are  determined  using  the  compatibility  equations,  [2.16],  and  [2.17]  and  the 
orifice  flow  equation,  [3.7]. 

Interface  3  (Feed  Line-Discharge  Line).  This  is  a  complex  boundary 
condition  which  models  the  combined  action  of  the  pump  and  its  inlet 


« 
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Table  4.2 


Monopropellant  System  Parameters 


Combustor  Properties 


Dead  Time,  Tc 

3 

msec 

Residence  Time, 

8  g  1 

msec 

Combustor 

Parameter 

370000 

1/ft2 

Flow  System  Properties 

Density,  p 

53 

lbm/ft3 

Steady  State 

Flow  Rate 

4 

ft'Vsec 

Injector  Head 
Drop,  AHinj 

40 

*/. 

Injector  Const., 

Binj  2- 8 

ft/sec 

Pump  Pressure 

Input  Motion  Amplitudes 

Rise  Coefficients 

Pump 

0.8 

ft/sec 

A' 

1208 

ft 

Tank 

0.32 

ft/sec 

B* 

668 

sec/ft2 

C* 

-109 

sec2/ft5 

- 

Sonic  Velocity 

Lenath  Area 

Friction  Factor 

(ft/sec ) 

<ft)  (ft2) 

Tank 

1000 

2  20 

0 

Feed  Line 

2000 

40  0.2 

0.032 

Discharge  Line 

3000 

6  0.05 

0.042 

Computer  Requirements 

Since  only  the  values  corresponding  to  the  current  time  need  be  stored 
(plus  a  few  old  values  and  a  small  queue),  the  implementation  of  the  routines 
developed  in  this  investigation  will  run  on  a  small  system  such  as  a  micro 
(desk  top)  computer.  In  fact,  all  of  the  programs  developed  in  this 
investigation  were  run  on  a  micro  computer  (the  Apple  Macintosh**). 

Single  Pine  Model.  This  particular  routine  requires  only  about  10  K 
(10,000  bytes  of  storage)  of  memory  for  its  source  code  and  about  2  K  for 
variables.  Typical  run  times  are  on  the  order  of  one  minute.  A  small 
system  is  well  suited  to  this  type  of  simple  model. 

Monoprooellewt  System.  This  routine  requires  much  more  storage 
(15  K  source  code  and  8  K  variables)  and  run  times  (for  transient  responses) 
are  about  two  hours.  This  type  of  performance  is  probably  adequate  for 
transient  resonses  (which  need  run  for  only  a  few  cycles  of  the  disturbance) 
but  would  not  be  adequate  for  steady  state  responses  involving  multiple  runs 
of  much  greater  duration. 

Bipropellant  System.  This  routine  requires  more  than  twice  the 
storage  of  the  monopropellant  system.  Transient  response  run  times  are 
about  six  hours  and  any  useful  par*.;»etric  study  should  be  carried  out  on 
either  a  'fast'  mini  computer  or  a  mainframe.  For  comparison,  this  program 
was  partially  implemented  on  the  ASD  CYBER  (CDC  6400)  main  frame 
computer.  Run  times  were  a  very  modest  200  CPU’s  (just  a  few  seconds  of 
real  time).  Future  users  of  these  routines  would  do  well  to  transport  them 
to  a  main  frame  computer  (see  Appendix  A). 
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V.  Results 


In  this  section,  results  from  the  single  pipe  and  monopropellant  system 
models  described  in  Section  IV  are  analyzed  and  compared  with  published 
results.  The  effects  of'various  parameter  choices  on  the  relative  POGO 
stability  of  the  monopropellant  and  bipropellant  models  developed  in  this 
investigation  will  also  be  examined. 


Single  Pine  Systems 

Fig.  5.1  is  a  plot  of  data  generated  by  the  single  pipe  program  described 
in  Section  IV.  Fig.  5.2  was  published  by  Wylie  and  Streeter14139  who 
originally  modeled  this  system.  Although  they  are  in  dissimilar  units,  a 
point  by  point  conversion  of  values  given  by  these  plots  shows  an  excellent 
agreement  between  the  results  generated  by  the  routine  developed  in  this 
study  and  Wylie  and  Streeter's  published  results.  A  number  of  other 
observations  follows  (note  that  pressure  and  velocity  indices  refer  to  Fig. 
4. IB). 

Valve  Flow  Velocity  (V1  D).  The  valve  flow  velocity  falls  off  starting 

at  time  zero  and  is  completely  stopped  in  2. 1  seconds,  according  to  the  valve 
flow  constant  described  by  {4.13  (although  the  increasing  valve  inlet 
pressure  modifies  the  shape  of  this  curve). 

Valve  Inlet  Pressure  (P,  y).  At  time  zero,  the  fluid  system  is  in 

steady  state  equilibrium  with  a  downstream  valve  inlet  pressure  of  1.41  MPa 
(143  feet  of  head  in  Fig.  5.2).  As  the  valve  closes,  a  large  pressure  rise 
(due  to  the  water-hammer  effect)  may  be  observed  in  both  plots.  After  its 
peak,  the  pressure  oscillates  with  a  period  of  about  2  seconds  which  is  the 
natural  period  of  the  system  (4  times  the  transit  time  of  the  pressure 
pulse). 
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Reservoir  Flow  Velocity  (VQ  y).  The  flow  velocity  at  time  zero  is  at 

the  steady  state  condition  of  2. 45  m/sec.  Due  to  the  finite  transfer  time  of 
the  pressure  pulse,  there  is  approximately  a  half-second  delay  in  the 
response  at  the  upstream  end  of  the  pipe.  The  flow  begins  to  fall  off, 
eventually  reversing  and  oscillating  about  the  zero  value  (at  the  system’s 
natural  frequency). 

Other  Single  Pine  Systems.  Since  the  object  of  this  thesis  is  to 
develop  and  document  routines  for  the  analysis  of  POGO  instability  in  liquid 
rocket  systems,  several  other  boundary  conditions  were  substituted  for  the 
simple  valve  and  constant  pressure  reservoir  of  Wylie  and  Streeter’s  model. 
POGO  instability  occurs  during  the  normal  full  throttle  operation  of  a  liquid 
rocket and  therefore  the  valve  transient  represented  by  [4.1]  was 
deleted  in  favor  of  a  constant  valve  flow  constant.  This  flow  constant  was  to 

be  set  to  T0,  the  initial  value  used  by  Wylie  and  Streeter  but  was 

unintentionally  set  to  500- 10“6  (m3/Kg),/2  which  is  25%  of  TQ. 

POGO  instability  is  caused  by  pressure  oscillations  (due  to  normal  thrust 
variations)  In  the  propellant  feed  system  of  a  liquid  rocket,  which  give  rise 
to  propellant  injection  variations,  thereby  causing  greater  thrust 
variations'-2.  In  an  attempt  to  simulate  this  phenomenon  in  the  single  pipe 

model,  the  upstream  pressure  of  the  pipe,  PQ  y  is  given  a  sinusoidally 

oscillating  variation  (of  frequency  0.5  Hz,  the  systems  resonant  frequency) 
from  its  mean  value.  This  variation  is  described  by  [4.41.  The  mean  value  is 

still  taken  as  PQ,  the  same  value  as  in  Wylie  and  Streeter's  model. 

A  hypothetical  accumulator  is  also  added  to  the  valve  inlet  to  simulate 
the  cavitation  compliance  at  the  injector  of  a  liquid  rocket  (see  Sections  III 
and  IV  for  details).  This  compliance  is  modeled  with  the  isothermal 
expansion  result  [3.15]  with  the  initial  bubble  pressure  taken  as  the  valve 
inlet  steady  state  pressure  and  the  initial  bubble  volumes  taken  as  0 
(effectively  zero  compliance),  0.01,  0.1  and  1  m3.  These  volumes  are 
representative  of  initial  volumes  used  in  the  monopropellant  system  of 


Dorsch,  et  al. 1 

Results.  This  single  pipe  compliant  system  is  at  the  steady  state 
condition  (flow  velocity  of  480  mm/sec  and  valve  pressure,  P1  y  of  1.45 
MPa)  when  the  0.5  Hz  pressure  oscillations  begin.  Figs.  5.3  and  5.4  are 
plots  of  the  pressure  and  velocity  at  the  valve  inlets  and  outlets  (  P,  u  and 

V,  D  ,  respectively)  of  systems  with  initial  cavitation  bubbles  as  noted. 

Increasing  amounts  of  valve  cavitation  compliance  (increasing  initial  bubble 
volume)  lead  to  decreased  pressure  and  velocity  perturbations.  Also,  the 
zero  compliance  responses  are  still  building  after  5  sec,  while  the  compliant 
system  responses  are  bounded.  This  is  the  main  purpose  of  the  compliant 
devices  (i.e.  accumulators  and  ‘fix‘^:3  devices)  introduced  in  the  pump 
inlets  of  most  liquid  rockets. 

Fig.  5.5  plots  steady  state  valve  inlet  pressure  disturbance 
amplitudes  (peak  to  peak)  for  initial  valve  inlet  cavitation  bubble  volumes  of 
0  m3  (zero  compliance)  and  1  m3.  For  this  plot,  the  inlet  pressure 

oscillation  frequency,  u  is  varied  from  0.1  Hz  to  3  Hz  (0.63  to  18.9 
rad/sec).  Note  the  effect  of  the  valve  compliance  (compliance  increases 
with  initial  accumulator  volume)  is  to  shift  the  resonance  phenomena  down  in 
frequency  and  to  greatly  decrease  its  overall  amplitude. 


The  monopropellant  system  model  as  detailed  in  Section  IV  was  run  with 
the  same  data  as  the  system  described  by  Dorsch,  et  al. '  in  order  to  verify 
the  routines  developed  in  this  investigation.  A  number  of  runs  were 
compared  with  those  already  published;  three  comparisons  will  be  presented 
here.  Before  any  meaningful  comparisons  can  be  made,  however, 
performance  criteria  must  be  established.  The  two  factors  most  important  in 
POGO  instability  are  pump  inlet  and  combustion  chamber  conditions.  Note  that 
the  indices  of  the  pressures  and  velocities  refer  to  Fig.  4.3B. 


Figure  5.3 

Single  pipe  valve  pressures  for  various  compliances . 


Figure  5.4 

Single  pipe  valve  flow  velocities  for  various  compliances. 
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Figure  5.5 

Frequency  responses  for  zero  and  larye  compliances. 


Pump  Inlet  (P3  y  and  V3  y).  The  performance  of  the  pump  is  directly 

related  to  the  pump  inlet  pressure'* '®( see  Table  4.2).  The  fluctuation  of  the 
pump  inlet  pressure  due  to  the  oscillatory  pump  and  tank  motions  will 
therefore  be  a  primary  indicator  of  POGO  instability  (pump  inlet  flow 
variations  will  be  of  similar  importance).  If  the  inlet  pressure  minima  drive 
the  pump  into  cavitation  (surge),  the  resulting  pump  output  pressure  will 
greatly  magnify  these  minima,  since  the  pump  pressure  rise  is  very  much 
greater  in  magnitude  (and  therefor  its  variation  will  be  much  greater)  than 
the  pump  inlet  pressure.  These  pressure  oscillations  will  lead  to  injector 
flow  velocity  variations  which  are  in  turn  fed  back  into  the  system  as 

variations  in  the  vehicle  accelerate, 

Combustor  (P«  D  and  V5  D).  As  previously  stated,  pump  output 

pressure  oscillations  give  rise  to  injector  flow  variations  which  are  fed  back 
into  a  real  system  (the  monopropellant  model  employs  a  constant  structural 
accleration  and  thus  ignores  this  feed  back)  in  the  form  of  structural 
acceleration  variations.  The  presence  of  large  combustor  pressure  and  flow 
oscillations  (driven  by  pump  output  pressure  oscillations)  indicates  POGO 
instability'1^. 

The  pump  inlet  and  combustor  states  were  determined  for  several  pump 
inlet  compliance  values.  In  each  case  examined,  the  agreement  between  the 
routines  developed  in  this  study  and  the  analysis  of  Dorsch,  et  al.  was 
striking. 

Zero  Compliance.  As  a  basis  upon  which  to  compare  later  runs, 
Dorsch,  et  al.'  examined  the  response  of  the  monopropellant  system  with  no 
accumulators  to  input  sinusoidal  pump  and  tank  motions  with  amplitudes  of 
0.8  ft/sec  and  0.32  ft/sec,  respectively.  The  excitation  frequency  was 
chosen  as  10  Hz. 

Results.  Fig.  5.6  plots  the  pump  suction  pressure  (P3  y)  transient 

response  of  the  monopropellant  system  as  modeled  using  the  routines 
detailed  in  this  investigation.  Fig.  5.7  plots  the  same  response  as  published 
by  Dorsch,  et  al.':'^.  Although  the  scales  and  units  of  the  two  plots  are 
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Figure  5.6 

Program  results  for  the  monopropellant  system  pump  inlet  pressure  (zero  compliance). 


Figure  5.7 

Pump  inlet  pressure  results  from  Dorsch,  et  al.1  (zero  compliance). 


dissimilar,  a  conversion  of  the  units  of  Fig.  5.7  to  those  of  Fig.  5.6  showed  a 
quantitative  agreement  (to  within  approximately  2X).  Note  also  the  excellent 
agreement  of  the  wave  forms,  both  in  the  initial  transient  and  in  the  final 
steady  state  amplitudes.  Only  a  non-linear  analysis  (i.e.  one  based  on  the 
Navier-Stokes  equations,  [2. 1]  and  [2.2])  could  have  reproduced  these  large 
amplitude  responses  since  the  pressure  amplitudes  are  far  to  large  as  to  be 
consistent  with  the  assumptions  used  in  a  linear  network  analysis* *4. 

Low  Compliance  Pumo  Inlet.  For  this  case,  the  pump  inlet  compliance 
was  modeled  with  the  isothermal  expansion  result  of  [3.15],  the  initial 
pressure  taken  as  1.47  MPa  (150  feet  of  ^0)  and  the  initial  volume  as  113 

cc.  The  applied  excitation  was  the  same  as  for  the  zero  compliance  case. 

Results.  Fig.  5.8  depicts  the  results  from  the  program  described  in 
this  Investigation,  while  Fig.  5.9  shows  published  results*524  for  the 

transient  responses  of  the  pump  inlet  velocity  (V3  y)  and  flow,  respectively. 

There  is  a  marked  similarity  between  these  very  non-linear  plots  (again, 
although  the  units  are  dissimilar,  they  correspond  very  closely  after 
conversion).  Apparently,  the  introduction  of  the  compliance  at  the  pump 
inlet  leads  to  very  complex  behavior  when  coupled  with  the  relative  pump 
motion  and  combustion  chamber  dynamics.  This  conclusion  is  supported  by 
the  fact  that  for  the  zero  compliance  case  these  same  plots  (not  shown)  were 
quite  sinusoidal.  Dorsch's  analysis  as  well  as  the  routines  developed  in  this 
investigation  provide  both  the  complex  transient  response  and  the  amplitude 
and  wave  form  of  the  steady  state  oscillations.  This  is  another  advantage  of 
finite  difference  based  methods. 

Lamp  Compliance  Pump  Inlet.  This  case  is  similar  to  the  previous 
low  compliance  case,  except  that  the  initial  bubble  volume  is  taken  as  283  cc, 
and  the  excitation  frequency  as  6  Hz.  The  effect  of  the  greater  inlet 
compliance  is  to  lessen  the  coupling  between  pump  inflow  and  outflow. 

Results.  Fig.  5. 10  is  a  plot  of  the  combustion  chamber  pressure 

(P5  0)  as  determined  by  the  routine  developed  in  this  investigation,  while 

Fig.  5.11  depicts  the  analogous  published  result***2^.  After  conversion  to 
similar  units,  the  plots  were  found  to  agree  within  approximately  ZX.  Notice 
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Figure  5.8 

Program  results  for  the  monopropellant  system  pump  inlet  velocity  (low  compliance). 


0 

.0* 

.» 

.12 

.1* 

.a 

.2* 

.21  .32 

.36 

.40 

Time,  sac 

Figure  5.9 

Pump  inlet  velocity  results  from  Dorsch,  et  al.1  (low  compliance). 
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Chanter  pressure  haa<  ft 


Figure  5.10 

Program  results  for  the  monopropellant  system  combustor  pressure  (large  compliance). 


Figure  5. 1 1 

Combustor  pressure  results  from  Dorsch,  et  al.1  (large  compliance). 
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the  relatively  sharp  nature  of  the  pressure  peaks.  This  is  due  to  the  pump 
characteristics.  As  the  inlet  pressure  drops,  the  pump  begins  to  cavitate 
producing  a  rapidly  decreasing  pressure  rise,  causing  the  pump  output 
pressure  to  exhibit  an  abrupt  pressure  minimum.  When  the  inlet  pressure  is 
near  its  peak,  the  pump  inlet  pressure  response  is  relatively  flat  leading  to  a 
more  rounded  pressure  maxima  in  the  pump's  output  pressure.  These  effects 
are  propagated  to  the  combustion  chamber  giving  rise  to  Figs.  5. 10  and  5.11. 


BiorooeMant  Liquid  Rocket 

As  a  guide  to  the  application  of  the  software  developed  in  this  thesis  by 

m 

future  investigators  of  more  complex  liquid  rocket  models,  a  typical  liquid 
rocket  system  was  examined.  The  closed  loop  (i.e.,  a  system  employing 
propulsion  feed  back)  model  was  described  in  Section  IV  while  its  actual 
implementation  in  terms  of  digital  computer  code  is  described  in  Appendix  A. 
In  this  section,  the  results  obtained  from  this  code  will  be  discussed.  The 
following  plots  are  important  in  determining  the  relative  P060  stability  of  a 
system  (note  that  the  indices  of  the  pressures  and  velocities  refer  to  Fig. 
4.6B). 

Pump  Inlet  Pressure  and  Velocity  (P,  u  and  V-.  u).  Just  8S  for  the 

monopropellant  system,  the  pump  inlet  pressure  variations  are  of  paramount 
concern  since  they  will  be  magnified  by  the  pump  inlet  pressure  response 
(see  Fig.  4.7).  The  pump  inlet  velocity  is  an  indicator  of  the  magnitude  of 
the  feed  line  now  variation,  as  well  as  the  relative  pump  structural  motion. 
The  presence  of  large  (i.e.  of  much  greater  magnitude  than  the  driving 
thrust  oscillations)  pump  inlet  pressure  and  velocity  oscillations  is 
indicative  of  a  P060  instability 1 : 1  ®. 

Combustor  Pressure  iM  Velocity  (P5  D  and  V5  0).  Since  the 

bipropellant  system  thrust  variations  will  be  fed  back  as  structural 
acceleration  variations,  the  injection  velocity  oscillations  are  of  great 
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importance.  The  combustor  pressure  is  also  examined  since  large  chamber 
pressure  variations  would  be  consistent  with  large  injection  oscillations  and 
since  other  modes  of  instability  can  be  excited  by  excessive  chamber 
pressure  variations ,:22. 

Component  Displacements  (XT  and  Xp).  During  a  POGO  event,  the 

tanks  and  pumps  will  oscillate  with  large  displacement  amplitudes3* ,0 
(presumably  at  or  near  their  natural  frequencies,  10  Hz  in  this  case).  These 
displacements  may  therefore  be  used  to  determine  the  presence  of  POGO 
instability. 

Both  the  transient  and  steady  state  response  (discussed  at  the  end  of  this 
section)  were  determined  for  three  different  configurations  with  varying 
degrees  of  pump  inlet  and  injector  dome  compliance. 

Svatcm  with  Zero  Compliance.  As  a  basis  for  comparison,  the 
compliances  of  the  pump  inlet  and  injector  dome  cavitation  bubbles  were 
taken  as  zero  (accumulator  initial  volumes  equal  to  zero).  This  implies  that 
the  pump  and  injector  outflows  equal  their  respective  inflows  (corrected  for 
relative  motion).  The  system  was  excited  with  a  sinusoidal  thrust  variation 
of  50,000  N  amplitude  at  a  frequency  of  12  Hz.  This  magnitude  represents  a 
modest  10%  thrust  variation  (as  defined  in  Appendix  B,  the  nominal  thrust  is 
500,000  N).  The  12  Hz  response  is  plotted  since  the  system  exhibited  a 
resonance  near  this  frequency  (see  steady  state  response). 


Pump  Inlet  Pressures  ( P,  y  F  and  P3  y  Qh).  Fig.  5.12  depicts  the 

transient  response  of  the  fuel  and  oxidizer  pump  inlet  pressures  to  the 
sinusoidal  thrust  excitation.  The  fuel  pump  inlet  exhibits  an  essentially 
sinusoidal  response,  while  the  much  less  reponsive  oxidizer  system  shows  a 
very  non-linear  response  due  to  its  lower  natural  frequency.  This  lower 
natural  frequency  occurs  because  the  oxidizer  lines  are  more  compliant  than 
the  fuel  lines  (see  Appendix  B),  the  larger  compliance  causing  a  lower  sonic 
velocity  (and  therefore  natural  frequency14)  as  defined  by  [4.7J.  Note  also 
that  both  the  fuel  and  oxidizer  pump  inlet  pressures  are  in  the  stable  (i.e,. 
flat)  region  of  Fig.  4.8.  This  indicates  the  relative  stability  of  the  system  to 
the  thrust  variation.  This  would  not  necessarily  be  true  if  different  pump 
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Single  Pine  Models.  The  results  from  a  single  pipe  model  utilizing  the 
/ 

omputer  routines  developed  in  this  investigation  compared  favorably  with 
>e  published  results  of  Wylie  and  Streeter14.  It  was  also  shown  that  valve 
ilet  compliance  (modeled  using  [3.15])  attenuated  the  effects  of  inlet 
ressure  disturbances. 

Honopropellant  System.  A  more  complex  system  first  modeled  by 
orsch,  et  al.1  was  modeled  using  the  software  developed  in  this  thesis.  A 
ery  favorable  comparison  with  results  published  by  Dorsch  verified  that  the 
outines  developed  in  this  study  produced  reliable  information. 

Bipropellawt  System.  A  typical  bipropellant  system  was  designed 
sing  a  simple  steady  state  analysis  (see  Appendix  B)  and  the  resulting  model 
/as  analyzed  using  the  routines  developed  in  this  investigation.  The 
ombination  of  a  small  excitation  amplitude  and  conservative  pump  operating 
haracteristics  yielded  results  that  contra-indicated  P060  instability, 
urthermore,  a  very  limited  parametric  study  of  the  effects  of  pump  inlet 
md  injector  dome  compliances  showed  this  system  to  be  very  insensitive  to 
hese  values  due  to  the  small  input  thrust  variation  and  stable  pump 
haracteristics. 


Figure  5.26 

Chamber  pressure  frequency  response  (zero  compliance) 
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System  with  Intermediate  Compliance.  This  system  produced  pump 
inlet  pressure  responses  (not  shown)  identical  to  those  for  the  zero  and 
large  compliance  cases. 

Combustor  Response.  Fig.  5.26  plots  the  combustor  pressure 
response  amplitudes  for  the  zero  compliance  case.  The  excitation  was  the 
same  sinusoidal  50,000  N  thrust  variation  used  throughout  this  section.  The 
curves  for  the  large  and  intermediate  compliance  cases  were  identical  to  the 
zero  compliance  case  and  are  not  shown  here.  The  zero  compliance  curve 
indicates  amplitude  (difference  between  maxima  and  minima)  minima  at 
approximately  11  Hz  and  14  Hz  with  larger  amplitudes  at  both  higher  and 
lower  frequencies.  Note  that  ail  amplitudes  are  very  much  smaller  than  the 
mean  chamber  pressure  (much  less  than  IX  of  mean  chamber  pressure). 
This  magnitude  is  considerably  smaller  than  that  for  the  pump  inlet  pressure 
(about  10X  of  mean  pressure),  giving  another  indication  of  the  stability  of  the 
pump  output  pressures. 


Figure  5.24 

Fuel  pump  inlet  pressure  frequency  response  (large  compliance). 


Figure  5.27 

Ox.  pump  inlet  pressure  frequency  response  (large  compliance) 


Figure  5.22 

Fuel  pump  inlet  pressure  frequency  response  (zero  compliance). 


Figure  5.23 

Ox  pump  inlet  pressure  frequency  response  (zero  compliance) 
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and  velocity  responses  do  not  necessarily  occur  at  component  boundaries.  If 
these  values  are  important  to  a  particular  analysis,  responses  at  points 
interior  to  a  fluid  lines  could  easily  be  obtained. 

Pump  Inlet  Response  (P3  j).  As  was  the  case  for  the  transient 

response,  the  steady  state  pump  inlet  pressure  response  is  of  great 
importance  since  low  pressures  can  cause  pump  inlet  cavitation  and  large 
pump  output  pressure  disturbances.  At  a  given  frequency,  the  presence  of 
large  pump  inlet  pressure  steady  state  amplitudes  resulting  in  pressure 
minima  lower  than  the  stable  operating  point  of  the  pump  (Fig.  4.8),  would 
indicate  a  POGO  instability' 

System  with  Zero  Compliance.  Figs.  5.22  and  5.23  depict  the  fuel 
and  oxidizer  pump  inlet  pressure  response  of  the  zero  compliance  system. 
On  each  plot,  the  upper  curve  represents  the  pressure  maxima  (i.e.,  the 
maximum  steady  state  pressure  exhibited  at  the  given  frequency),  while  the 
lower  curve  reflects  the  pressure  minima  (i.e.,  the  minimum  steady  state 
pressure).  The  dashed  line  between  the  curves  is  the  unperturbed  (i.e,. 
initial)  pressure.  Note  that  the  very  non-linear  nature  of  the  plot  is  due  in 
part  to  the  small  number  of  data  points  taken  ( 15  per  plot).  The  fuel  pump 
traces  have  peaks  at  8  and  12  Hz,  and  approach  the  mean  value  (nominal 
pump  inlet  pressure  of  180  KPa)  at  high  and  low  frequencies.  The  oxidizer 
pump  inlet  response  is  more  complex  and  seems  to  indicate  peaks  near  6,  13, 
and  22  Hz,  with  decreasing  response  (difference  between  maxima  and 
minima)  at  higher  and  lower  frequencies.  In  neither  case  do  the  pressure 
minima  approach  the  cavitation  points  in  Fig.  4.7,  thus  it  can  be  concluded 
that  the  system  exhibits  no  POGO  instability  over  the  range  of  frequencies 
examined. 

System  with  Large  Compliance.  Figs.  5.24  and  5.25  are  the 
analogous  plots  for  the  system  with  283  cc  and  50  cc  initial  pump  inlet  and 
injector  dome  accumulator  volumes  respectively.  These  plots  are  virtually 
identical  to  the  zero  compliance  plots  (Figs.  5.24  and  5.25),  indicating  the 
insensitivity  of  the  pump  inlet  pressures  to  inlet  compliance  values  at  this 
level  of  input  thrust  oscillation. 


Figure  5.20 

Bipropellant  system  pump  inlet  pressures  (intermediate  compliance). 


msec 

Figure  5.21 

Bipropellant  system  combustor  pressure  (intermediate  compliance). 
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System  with  intermediate  Compliance.  For  this  case,  the  initial 
bubble  volumes  for  the  fuel  and  oxidizer  inlets  were  lowered  to  140  cc,  while 
the  fuel  and  oxidizer  injector  dome  values  were  raised  to  250  cc.  Although 
the  resulting  data  reflects  both  of  these  changes  (making  it  difficult  to  judge 
their  effects  separately),  the  system  has  proven  to  be  relatively  insensitive 
to  the  localized  compliance  values.  This  model  produced  responses  identical 
to  the  zero  and  large  compliance  cases.  Two  plots  are  given  for  comparison. 

Pump’lnlet  Pressures.  Fig.  5.20  is  the  graph  of  the  fuel  and  oxidizer 
pump  inlet  velocites  for  the  intermediate  compliance  case.  They  appear 
almost  identical  to  the  previous  two  cases,  indicating  the  relative 
insensitivity  of  the  inlet  pressures  to  the  compliance  values  (at  this  level  of 
Input  disturbance). 

Combustor  Pressure.  The  combustor  pressure  is  given  in  Fig.  5.21. 
Again  the  pressure  fluctuation,  though  very  complex  in  wave  form,  is  very 
small.  This  particular  system  is  indeed  quite  stable. 

State  Response 

A  principal  advantage  of  the  method  of  characteristics  based  approach 
used  In  this  analysis  is  the  ability  to  obtain  both  the  transient  snd  the  steady 
state  responses.  The  steady  state  response  is  simply  the  final  waveform  and 
amplitude  of  the  transient  response  plot  after  the  transient  has  died  out. 
This  principle  can  be  applied  at  all  points  in  the  system  for  arbitrary 
frequencies,  resulting  in  the  complete  steady  state  description  of  the 
system's  response  to  any  input.  However,  due  to  the  computational 
limitations  of  this  analysis,  only  a  small  number  of  plots  (pump  inlet  and 
combustor  pressures)  will  be  examined.  Furthermore,  only  a  small  number 
of  points  (each  corresponding  to  the  response  at  a  discreet  frequency)  will 
be  taken.  It  should  be  noted  that  this  is  not  an  inherent  limitation  of  the 
methods  presented  in  this  thesis  and  a  more  complete  analysis  of  the 
pressure  and  velocity  steady  state  responses  at  many  more  points  in  the 
system  would  be  most  useful.  Note  that  the  largest  steady  state  pressure 


Figure  5.10 

Bipropellant  system  pump  inlet  pressures  (large  compliance). 


msec 


Figure  5. 19 

Bipropellant  system  combustor  pressure  (large  compliance). 
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is  about  0.5  mm  which  represents  a  10%  variation  from  the  mean  value  (5 
mm).  Since  the  original  thrust  variation  was  10%,  a  stable  system  is 
indicated^  *3»4. 

Structure  Acceleration.  Fig.  3.17  gives  the  structural  acceleration 
of  the  zero  compliance  system.  Note  that  the  structure  includes  all  of  the 
system  components  except  the  tanks  and  pumps  (see  Fig. 4. 7).  A  large 
start-up  transient  of  approximately  80  msec  duration  is  present.  This  is  due 
to  the  small  transient  response  of  the  component  displacements,  magnified 
by  the  fact  that  the  structure  is  very  light  compared  to  the  tanks.  As  in  the 
case  of  the  component  displacements,  the  relative  magnitude  of  the  steady 
state  structural  acceleration  is  the  same  as  the  original  disturbance  (10%) 
indicating  the  system's  stability. 

System  with  Large  Compliance.  This  system  is  the  same  as  the  zero 
compliance  system  except  that  the  pump  inlet  and  injector  dome  compliances 
were  modeled  with  the  isothermal  expansion  model  (see  Section  II.).  The 
initial  bubble  pressures  were  chosen  as  the  steady  state  (nominal)  pressures 
(Table  4.5).  The  initial  fuel  and  oxidizer  inlet  bubble  volumes  were  each  Z83 
cc,  while  the  initial  fuel  and  oxidizer  injector  dome  bubble  volumes  were  50 
cc.  This  model  produced  results  that  were  almost  identical  with  the  zero 
compliance  case  indicating  that  the  initial  thrust  disturbance  was  too  small  to 
cause  any  significant  pump  ouput  variations,  and  thus,  no  feed  back  in  the 
form  of  flow  velocity  variations  reached  the  combustion  chamber.  Two 
typical  plots  are  given  below. 

Pump  Inlet  Pressures  (P3  Ut  f  and  ^3  u  Ox^*  5* ,8  detai,s  the 

time  history  of  the  fuel  and  oxidizer  pump  inlet  pressures.  There  is  a  great 
similarity  to  the  zero  compliance  case  (Fig.  5.14),  although  Fig.  5.18  is 
slightly  smoother  In  comparison.  The  localized  pump  inlet  compliance  seems 
to  have  little  effect  on  this  very  stable  system. 

Combustion  Chamber  Pressure  (P5  D  F  and  (P5  D  0k).  The 

combustor  pressure  response  is  given  in  Fig.  5.19.  Again  this  curve  is 
almost  identical  to  the  zero  compliance  case,  the  resulting  oscillation  being 
very  small.  This  indicates  the  lack  of  POGO  instability. 


2  uP 3 


Figure  5.14 

Bipropellant  system  combustor  pressure  (zero  compliance) 


Figure  5.15 

Bipropellant  system  injector  velocities  (zero  compliance). 
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Figure  5. 1 3 

Bipropellant  system  pump  inlet  velocities  (zero  compliance). 
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characteristics  had  been  chosen  in  Appendix  B.  The  importance  of  accurate 
pump  characteristics  cannot  be  over  stated. 

Pump  Inlet  Velocities  (V3  u  F  and  V3  y  Fig.  5. 13  is  a  plot  of 

the  pump  inlet  velocities  for  the  zero  compliance  system.  As  in  the  case  of 
the  inlet  pressures,  the  fuel  pump  inlet  velocity  trace  is  much  more 
sinusoidal  than  the  very  non-linear  oxidizer  pump  Inlet  velocity  trace 
(another  manifestation  of  the  oxidizer  systems  lower  sonic  velocities).  The 
fuel  pump  velocity  amplitude  is  approximately  3%  of  the  nominal  value  while 
the  corresponding  oxidizer  amplitude  is  7%  of  its  nominal  value.  These  very 
modest  variations  further  indicate  a  stable  system. 

Combustion  Chamber  Pressure  ( P5  D  p  and  ( P5  D  ^ ) .  The  chamber 

pressure  transient  response  is  given  in  Fig.  5.M.  This  plot  is  much  more 
complex  and  non-linear  than  those  for  the  pump  inlets  due  to  the 
superposition  of  the  combustor  dynamics  on  the  pump  ouputs.  Note  also  that 
the  magnitude  of  the  chamber  pressure  oscillation  is  quite  small  compared  to 
the  mean  pressure  (another  indication  of  the  system's  stability). 

Injector  Outlet  Velocities  (V5  D  p  and  (V5  Q  The  fuel  and 

oxidizer  injector  outlet  velocities  for  the  zero  compliance  system  are  given 
in  Fig.  5.15.  While  both  are  quite  non-linear,  these  perturbations  are  a 
much  smaller  fraction  of  the  nominal  values  than  were  the  corresponding 
pump  inlet  velocities,  attesting  to  the  very  stable  operation  of  the  pump.  The 
pump  is  essentially  isolating  the  dishcharge  lines  and  combustor  from  the 
modest  pump  inlet  pressure  fluctuations. 

Component  Displacements  (XT  F).  For  simplicity,  the  pump  and  tank 

masses  and  structural  coefficients  were  selected  to  produce  dynamically 
identical  responses  with  resonances  at  10  Hz  (see  Table  4.4  and  Appendix  B). 

Fig.  5.16  plots  the  displacement  of  the  fuel  tank,  XT  p  in  response  to  the 

excitation  (the  oxidizer  tank  and  both  pump  displacements  have  responses 
identical  to  the  fuel  tank  responses).  The  transient  is  very  short  and  the 
steady  state  wave  form  is  nearly  triangular.  The  amplitude  of  the  response 
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VI.  Conclusions 


This  section  summarizes  the  findings  of  this  study  and  recommends  areas 
for  future  work.  The  accuracy  and  utility  of  the  routines  developed  in  this 
paper  in  modeling  POGO  instability  problems  of  arbitrary  complexity  has  been 
shown.  The  following  has  been  accomplished: 

1 )  A  general  tool  has  been  developed  for  the  study  of  liquid  rocket  POGO 
instability. 

2)  The  routines  developed  in  this  paper  have  been  verified  by  favorable 
comparison  with  published  data  from  several  sources. 

2)  A  systematic  reference  for  future  POGO  stability  analyses  has  been 
presented.  Core  routines  have  been  written  and  are  available  in  the 
Appendix. 

Recommendation 

The  model  developed  here,  while  complete,  could  be  extended  primarily 
with  the  addition  of  more  accurate  component  models.  Furthermore,  a 
complete  parametric  study  would  be  useful  in  actually  evaluating  the  effects 
of  various  components  on  the  overall  stability  of  a  given  system. 
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Appendix  A:  A  User's  Guide  to  PQGO  Analysis  Software 


This  appendix  Is  intended  as  a  guide  to  applying  the  routines  developed  In 
this  paper  to  arbitrary  liquid  rocket  systems.  Although  it  is  true  that  a 
completly  general  model  of  something  as  complex  and  specific  as  a  liquid 
rocket  system  Is  not  possible,  It  Is  possible  to  develop  core  routines  to  be 
used  in  the  modeling  and  simulation  of  arbitrary  systems.  The  routines 
documented  herein  should  be  viewed  as  an  aid  to  the  analyst  in  constructing  a 
complete  model,  not  as  such  a  model  in  themselves.  In  this  spirit,  various 
component  modeling  routines  will  be  presented  and  a  general  framework  for 
their  application  will  be  given. 


Cpmputational  Requirements 

The  routines  used  developed  in  this  investigation  and  detailed  in  this 
appendix  do  not  require  a  great  deal  of  computational  power.  This  is  due  to 
the  relatively  large  time  steps  possible  with  the  method  of  characteristics 
fluid  flow  solution.  Also,  since  only  one  set  of  conditions  must  be  stored  at 
any  one  time,  the  necessary  data  storage  In  arrays  and  queues  is  quite 
small.  Transient  solutions  can  be  obtained  conveniently  (with  four  or  five 
hours  of  computational  time)  on  even  a  micro  computer.  The  calculation  of 
steady  state  solutions  will  require  many  more  times  this  computations!  time 
and  mini  or  main  frame  computer  must  be  used. 

The  code  presented  here  is  written  in  the  MS-BASIC"'  language  (the  actual 
dialect  is  for  the  Apple  Macintosh"'  personal  computer).  This 
implementation  of  BASIC  is  sufficiently  general  such  that  it  will  run  with 
minor  or  no  modification  on  most  mini  and  micro  computers.  Many  main 
frame  systems  also  have  BASIC  interpreters,  although  some  modification  of 
the  routines,  especially  in  the  file  input-output  blocks,  will  likely  be 
necessary.  Note  that  this  version  of  BASIC  assumes  double  precision  (14 
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digits)  unless  otherwise  specified  by  the  following  rule:  variables  ending  in 
'X*,  T,  and'#'  are  of  integer,  single,  and  double  precision,  respectively. 


Single  Pipe  Systems 

AH  systems  will  consist  of  a  collection  of  single  pipe  elements.  The 
complexity  arises  from  the  combination  of  a  large  number  of  these  elements. 
A  simple  routine  to  evaluate  the  pressure  and  velocity  at  all  points  internal 
to  an  idealized  pipe  follows. 

Initialization  Block.  Before  the  routine  can  be  started,  the  following 
initialization  statements  must  be  executed: 


10  REH  Single  Pipe  Sgetea. 

40  GOSUB  200  ‘Initialization. 

41  GOSUB  225  'Calculate  Steady  State  Conditions. 

42  SVST 1 11E— 0  'Begin  Calculations. 

Body.  The  following  code  will  actually  apply  the  initial  boundary 
conditions,  calculate  the  state  at  the  internal  sections  of  the  pipe,  and  apply 
the  terminal  boundary  conditions. 

44  IF  I NT(SYST I IIE/DELTflTOUT) > I HT ( (SVST I OE-DELTRT)/DELTRTOUT) 

THEN  GOSUB  300  'Output  Results. 

45  SVST  I  IIE=SYST  I HE+DELTRT  '  I  ncreient  T  i  ae . 

46  IF  SYSTH1E>ENDTIt!E  THEN  END  'Stop  if  Coapleted. 

50  GOSUB  600  'Rpply  Initial  Boundary  Condition. 

60  FOR  IDOX-1  TO  NSECTIONSX-1  'Step  through  all  internal 
sections. 

70  GOSUB  400  'Get  Cp. 

73  GOSUB  500  'Get  Ca. 

75  OLOPRES-PRESSURE(IDOX):  OLDUEL=UELOCITV( IDOX)  'Save  old 
state  (for  next  section). 

80  PRESSURE (I BOX )»(CP+CH)/2  'Calculate  nee  pressure. 

90  UELOC I  TV ( I DOX ) * ( PRESSURE (ID0X)-C!i)/( RH0*S0UND )  'Calculate 

nee  velocity. 


101  NEXT  I DOS 

105  G0SU6  700  'Apply  Terminal  Boundary  Conditon. 

110  GOTO  44  'Repeat  until  finished. 

Initialization.  Before  a  program  can  run,  values  must  be  assigned  to 
the  necessary  variables.  The  data  read  in  will  be  defined  later  in  this 
appendix.  This  section  also  computes  the  number  of  reaches  necessary  to 
meet  the  Courant  condition  at  the  specified  time  step  size.  This  fixes  the 
dimension  of  the  storage  arrays  the  program  will  require. 

200  REff  Get  data  and  initialize  routine. 

205  READ  0ELTHT,  TEND,  DELTA TOUT,  RHO,  SOUND,  LENGTH,  DIRfl, 

RRER,  FRICT,  TRUO,  PBflSE,  RNP,  FREQ,  UGUES,  RCCEL,  RLPHR, 

TCUT,  TRUEXP 

210  HSECT I ONSX51 1 NT ( LENGTH/ (DELTRT*(SOUND+UGUES) ) )  'Use  Courant 

Condition  to  calculate  the  necessary  nuaber  of  reaches  for 
the  pipe. 

215  THETH=OELTHT/(LEHGTH/NSECT I ONSX) :  PS1=THETR*S0UND 
'Calculate  Grid  Parameters. 

220  Dill  UELOCITY(NSECTIONSX):  DM  PRESSURE(NSECTIONSX) 

'Dimension  state  arrays. 

222  RETURN 

Steady  State  Initialization.  The  finite  difference  routine  requires 
starting  values  on  which  to  base  its  calculations.  In  this  section,  the  steady 
state  response  is  calculatated  with  completely  conventional  methods.  Note 
that  the  steady  state  calculations  depend  upon  the  boundary  conditions.  This 
particular  section  assumes  an  initially  constant  pressure  upstream,  and  a 
valve  downstream  (orifice  with  decreasing  flow  coefficient). 

225  REN  Steady  State  Initialization. 

230  UEHD= ( ( TRU0*2* ( PBRSE+RHO*RCCEL*LENGTH*S I N(RLPHfl) ) )/ 

( 1 +(RH0*FR  I  CT*LENGTH*TRU0''2 ) / ( 2*0 1  fin) ) ) * .  5  'Calculate 

outlet  velocity. 

235  PEND»PBRSE+RHO*RCCEL*LENGTH*S I H(RLPHR)-(RHO*FR I CT*LENGTH* 
UEWT2)/  (2*0IRfl)  'Calculate  outlet  pressure. 

240  FOR  IDOX-O  TO  HSECT I ONSX 

250  PRESSURE( I DOX)“PBflSE+(PEND-PBflSE)* I OOX/NSECT I ONSX 


'Interpolate  interior  pressures. 

260  UEL0CITY(ID0X)»UEND  'Assign  constant  interior  velocity. 
270  NEXT  I  DOS 
275  RETURN 


Output.  This  section  will  of  course  vary  with  the  user's  wishes.  The 
following  routine  simply  prints  out  the  terminal  conditions  (in  scientific 
notation  to  5  significant  figures). 

300  REN  Output  Routine. 

310  PRINT  USING  "♦*.«*•* - _  SYSTIflE, 

PRESSURE(O) jUELOC I TY(0) , PRESSURE(NSECT I 0NSX) , UEL0C I TY(NSECT I 
QNSX) 

320  RETURN 


Determination  of  Cp  and  Cm.  This  section  calculates  the  values  of  Cp 
and  Cm  at  specified  sections  of  the  pipe. 

400  REN  Get  Cp. 

410  UELR»(UEL0C 1 TY( I D0X)-PS I *(UEL0C I TY( I 00X)-0LDUEL) ) 

/( 1 +THETR*(UEL0C I TY( I DOX) -OLDUEL) )  ' Co I cu I ote  Ur . 

420  PRESR»PRESSURE( I 00X)-(UELR*THETR+PS i )*(PRESSURE( I DOX) 

-0L0PRES)  ‘Calculate  Pr 

430  CP=PRESR+RH0*UELR*(S0UND+RCCEL*S I N ( ALPHA ) *DELTRT*S0UND/ 
UELR-S0UN0*FR I CT*DELTRT*RBS(UELR)/(2*D I RN) )  'Assign  Cp 

value. 

440  RETURN 

500  REN  Get  Ca. 

510  UELSs(UEL0C I TY( I 00X)-PS 1 *(UEL0C I TV( 1 D0X)-UEL0C I TV( I D0X+1 ) ) ) 

/  ( 1 -THETfl*(UEL0C I TV( I D0X)-UEL0C I TV( I D0X+1 ) ) )  'Calculate 

Us. 

520  PRESS=PRESSURE( I D0X)+(UELS*THETR-PS I )* 

(PRESSURE ( I D0X)-PRESSURE( I D0X+1 ))  'Calculate  Ps 

530  CN-PRESS-RH0*UELS*(S0UN0+RCCEL*S I N(RLPHfl)*DELTRT*S0UND/ 
UELS-S0UHD*FR I CT*DELTRT/(2*D I RN) )  'Assign  Ca  value. 

540  RETURN 


Boundary  Conations.  The  only  remaining  issue  is  the  application  of  the 
various  boundary  conditions  used  to  model  the  end  conditions  of  the  pipe.  One 


initial  condition  (known  pressure)  and  two  terminal  conditions  (valve  closure 
with  and  without  compliance)  will  be  presented. 

Initial  Condition.  The  initial  boundary  condition  for  this  system  is 
that  of  a  known  (oscillatory)  pressure.  The  following  code  matches  the 
known  upstream  pressure  to  the  downstream  conditions  (based  on  Cm). 

600  REff  Initial  Boundary  Condtion. 

610  I00X-0:  GOSUB  500  ‘Get  Cm. 

615  OLDUEL-UELOCITV(O):  OLOPRES-PRESSURE(O) 

620  PRESSURE ( 0 ) -PBRSE+RHP*S I N( 6 . 2832*FREQ*SVST I RE )  ‘Calculate 

upstreaa  pressure. 

630  UELOCITV(O)-(PRESSURE(0)-CH)/(RHO*S0UND)  ‘Hatch  velocity 
to  doenstreaa  conditions. 

640  RETURN 


Valve  Closure  with  Zero  Compliance.  The  exponential  valve  closure 
is  applied  to  the  outlet  flow  coefficient  and  the  conditions  are  matched. 

700  REN  Terminal  B.C. 

710  ID0X=NSECTI0NSX:  GOSUB  400  'Get  Cp. 

712  0LDPRES=PRESSURE( NSECT I 0NSX) :  0LDUEL=UEL0C I TV(NSECT I 0NSX) 

'Save  old  state  for  next  iteration. 

715  IF  SVST I tlE<TCUT  THEN  TflU-TfiUO* ( 1  -SVST I flE/TCUT ) ATRUEXP  ELSE 
TRU<*0  'Calculate  floe  coefficient. 

720  PRESSURE (NSECT I 0NSX)*( ( <RH0A2*S0UNDA2*TRUA2+4*CP) A . 5- 
RHO*SOUND*TRU)/2)A2  ‘Natch  pressure. 

730  UEL0C I  TV (NSECT 1 0NSX)»(CP-PRESSURE(NSECT 1 0NSX) )/(RH0*S0UND) 
'Calculate  velocity. 

740  RETURN 

Valve  Clossure  with  Compliance.  This  boundary  condition  routine 
allows  the  valve  inflow  and  outflow  to  differ  (using  the  compliance  relations 
developed  previously),  integrates  the  resulting  pressure  rate  of  change  (at 
averaged  conditions),  and  matches  the  resulting  conditions. 

700  REn  Terminal  B.C. 

710  I00X-NSECTI0HSX:  GOSUB  400  'Get  Cp. 

711  CPN I D"(0LDCP+CP)/2  ‘Calculate  average  Cp. 


78 


712  OLDPRES*PRESSURE(NSECTIONSX):  til DPRES=OLOPRES : 

OLOUELsUELOC I TV(HSECT 1 ONSX)  1  Guess  at  end  cond i t i ons . 

715  IF  SVSTIHE<TCUT  THEN  TflU=TRUQ*( 1 -SVST I nE/TQUT)ATflUEXP  ELSE 
TR(J3G  'Calculate  floe  coefficient. 

720  PRESSURE (HSECT 1 0NSX)=0LDPRES+DELTflT*(f1l DPRESi'2*flRER*( (CPU 1 0 
-H I DPRES)/  (RH0*S0UND)-TflU*(H 1 DPRES )A . 5)/(PH0T*NUH0T) ) 
'Integrate  pressure  rate  of  change. 

721  IF  RBS ( (PRESSURE ( HSECT I ONSX ) +0LDPRES ) /2-H I DPRES ) >T0L 1  THEN 
tl  I DPRES® ( PRESSURE ( HSECT I ONSX ) +0LDPRES ) /2 :  GOTO  720  'Check 
accuracy  of  guess,  repeat  if  insufficient. 

730  UELOC I TV ( HSECT 1 ONSX ) *TRU* ( PRESSURE ( HSECT I ONSX ) ) '".  5  'Hatch 

uekcity. 

740  RETURN 


Application  a  iynlcal  LiauM  BocKfil 

In  this  section  the  code  used  to  model  a  bipropellant  liquid  rocket  (see 
previous  sections)  will  be  presented  and  analyzed.  The  purpose  of  this 
section  is  to  facilitate  the  application  of  the  methods  and  software  developed 
in  this  paper  to  other  systems.  The  reader  will  note  the  great  similarity 
between  the  routines  for  this  complex  system  and  those  for  the  simple  pipe. 
The  essential  difference  is  that  most  simple  variables  from  the  last  program 
are  arrays  in  this  case.  A  detailed  list  of  all  program  variables,  both  those 
used  for  the  reading  of  data  and  those  internal  to  the  routine,  will  be 
presented  at  the  end  of  this  section. 

Dynamic  Response  Routine.  The  complete  liquid  rocket  model  consists 
of  two  routines:  a  dynamic  routine  employing  finite  difference  methods,  and  a 
steady-state  routine  to  initialize  the  finite  difference  model.  The  dynamic 
routine  simply  reads  the  initial  values  and  other  parameters  from  the  text 
(ASCII)  file  ’Steady  Results*.  The  program  text  follows: 

10  REH  BASIC  PR0GRRH  P0G0. 

15  C0HH0N  N0RTX,  DRTFILESO,  NDAT2X,  0RTFILE2$(),  PFILES, 

REH0TEX  'HI  loos  reaote  calling  froa  another  routine  (to 
gather  aany  steady  state  data  points  autoeatical ly). 


20  GOSUB  345  ‘DEFINE  VARIABLES  AND  RERO  STEAOV  STATE  RESULTS. 

25  GOSUB  170  'Get  user  input. 

30  GOSUB  1430  'Initialize  all  functions. 

35  REN  BEGIN  ROUTINE. 

40  IF  I NT(SVST I HE/DELTRTOUT) > I NT ( (SVSTI HE-DELTRT)/DELTRTOUT) 
THEN  GOSUB  1300  'Print  output  file. 

45  SVST I HEsSVST I NE+DELTRT :  IF  SVST I flE>TEHD  THEN  1375  'If 
completed  then  quit. 

50  GOSUB  945  'Establish  chaiber  properties. 

55  GOSUB  1200  'Establish  structural  response. 

60  FOR  I00X*FUELX  TO  0XIDI2ERX 

65  ON  (IDOX+1)  GOSUB  650,650  ‘Boundary  Conditions  Zero 
(ullage  conditions). 

70  FOR  DELNX=1  TO  NELNSX(IDOX) 

75  ELENENTX*QELNX 

80  GOSUB  630  'Get  PS I  and  THETR  (dependant  on  local  sonic 
uelocity). 

85  GOSUB  120  'Step  through  sections  (finite  difference 
routine). 

90  IF  ELENENTX=NELNSX( I DOX)  THEN  100  'If  completed,  goto 
teriinal  B.C. 

95  IF  I00X-FUELX  THEN  OH  ELE11ENTX  GOSUB  685,740,810,740  ELSE  ON 
ELEI1ENTX  GOSUB  685,740,810,740  'These  are  the  boundary 
condition  subrout i nes . 1 00  NEXT  DELHX  'Get  next  eleient. 

105  NEXT  I00X  'Goto  next  subsystem  (fuel  or  oxidizer)  or 
terminate  finite  difference  portion. 

110  GOSUB  1030  'Chaiber  Boundary  Conditions  (apply  teninal 
boundary  condition). 

115  GOTO  35  'Repeat  for  next  tiie  step. 


120  REN  Finite  difference  routine. 

125  FOR  I 0 1 X*NSECT I ONSX ( I DOX , ELENENT X ) + 1  TO 

NSECT I 0NSX( I DOX, ELENENTX+1 )-2  ' Increment  through  al I 

internal  sections  of  the  present  eleient. 

130  INTPRES-PRESSURE(IDIX):  1 HTUEL-UELOC 1 TV( I D1 X)  'Saoe 
intenediate  pressure  and  uelocity  for  next  section. 
135  GOSUB  580  'Get  CP. 

140  GOSUB  605  'Get  CN. 

145  PRESSURE( ID1X)*(CP+CH)/2 
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150  UELQCITV(  ID1X)«(PRESSURE(  ID1X)-Cn)/(S0UND(  1 DOX,  ELEI1EHTJE) 
♦RHO(IDOX)) 

155  OLOPRESC I DOX)= I NTPRES :  OLOUEL( I00X)=1NTUEL  'Pass  values. 
160  NEXT  I 01 X  'Goto  next  section. 

165  RETURN 

170  REN  Get  user's  input. 

175  IF  RENOTEX  THEN  OPEN  ’INSTRUCTIONS"  FOR  INPUT  RS  *99  'Check 
for  reaote  operation. 

180  IF  RENOTEX  THEN  INPUT  *99,  FREQ,  TSTRRT,  TEND,  DELTRTOUT, 
NORTX:  GOTO  230  'Read  reaote  operating  instructions  and 
branch. 

185  CLS:  CRLL  TEXTFONT(O):  CALL  TEXTS1ZE(18)  'Just  foraatting 
stuff. 

190  PRINT:  PRINT  TRB(11);"P0G0  RNRLVSIS  ROUTINE’ 

195  CALL  TEXTSIZEC12):  PRINT:  PRINT  TRB(5);"0riue  Frequency:  ’;: 
INPUT;  ”,FREQ 

200  PRINT:  PRINT  TRB(5); "Start  Tiae:  ’;:  INPUT;  ”, TSTRRT:  PRINT 
’  End  Tiae:  ’;:  INPUT  ”,  TEND 
205  PRINT:  PRINT  TRB(5); ’Output  Tiae  Increaent:  [  >“;DELTRT;’] 
INPUT  ”,  DELTRTOUT 

210  IF  DELTRTOUT<DELTRT  THEN  205  'Branch  if  choice  of  tiae 
increaent  is  unacceptable. 

215  PRINT:  PRINT  TRB(5);’This  sill 

generate"; INT((TEHD-SVSTINE)/DELTRT0UT)+1; "data  points.  OIC 
?  (V/N)’:  RES$*INPUT$(1)  'fll los  user  to  confira  before 
it's  too  late. 

220  IF  RESSO’V"  AND  RES$o’y’  THEN  170  'If  error,  try  again. 
225  PRINT:  PRINT  TRB(5);"Nuaber  of  data  series  to  generate 
(excluding  TINE  8  RCCEL):  ’;:  INPUT  ”,NDRTX 
230  Oin  DRTPX(NDRTX) :  DIN  DRTFILE$(NDRTX):  NDRT2X=1 :  DIN 
DRTF I LE2$(5) 

250  OPEN  *0",  *1 , "Tiaes’,32  'Open  and  initialize  ouput  files. 
255  OPEN  ’O’,  *2,"flccels’,32  'Initialize  structural  ouput 
f  i  les. 

260  OPEN  ’O’,  *3, ’Fuel  Tank  0isp.",10:  CLOSE  *3 
265  OPEN  ’O’,  *3,  "Fuel  Puap  Disp.’JO:  CLOSE  *3 
270  OPEN  "O’,  *3,  ’Ox.  Tank  Disp.’,10:  CLOSE  *3 
275  OPEN  ’O’,  *3,  ’Ox.  Puap  0isp.’,10:  CLOSE  *3 
280  FOR  IDOX-1  TO  NORTX  ' 


285  IF  RENOTEX  THEN  INPUT  *99,  SIDEX,  ELEHENTX,  SECX:  GOTO  300 
‘Get  reaote  particulars. 

287  REN  Get  inforaation  on  each  ouput  file  to  be  generated. 

290  PRINT:  PRINT  TRB(5);"SIDE:  (F/0)  ";:  RES$=INPUT$(1):  IF 
RES$="F"  OR  RES$="f"  THEN  SIDEX=FUELX  ELSE  SIDEX=0XIDI2ERX 
295  PRINT  RES$;"  ELENENT:  ";:  INPUT; “",ELENENTX:  PRINT  " 
SECTION:  INPUT;  "",SECX 

300  0RTPX( I 00X)=NSECT I ONSXCS I OEX, ELENENTX ) +SECX 
305  IF  RENOTEX  THEN  INPUT  *99,  DRTFILEX(IDOX):  GOTO  315 
310  PRINT  "  FILE:  ";:  INPUT  '  ”,DRTFILE$(IDOX) 

315  OPEN  "0",*2*ID0X+1,DRTFILE$(I00X)+\PRES":  CLOSE  *2*IDQX+1 
‘ Initial ize  ouput  files. 

320  OPEN  "0“ , *2* I DOX+2, DflTF I LE$( I 00X)+" . UEL " :  CLOSE  *2*ID0X+2 
325  NEXT  IDOX 

330  CLS  ’Clear  the  screen. 

335  IF  RENOTEX  THEN  CLOSE  *99 
340  RETURN 

345  REN  STEROV  STRTE  RESULTS  (input  the  steady  state  conditions 
as  detailed  in  the  file  "Steady  Results”. 

347  REN  Di aens ion  arrays. 

350  OIN  OLDCNI(I):  DIN  OLOCPI(I):  DIN  OLDUEL(I):  DIN  OLDPRES(J) 
DIN  0LDUELK1):  DIN  RH0(1):  OIN  ULLRGE(I):  DIN  NPUNPX(I): 
DIN  RRERINJECT(I):  DIN  TRUUI(I):  DIN  TRUU2(1):  DIN 
TRUINJECT(I):  DIN  PNOTI(I):  DIN  NUNOTI(I):  DIN  TDENDO(I): 
DIN  TUENDO(I) 

355  DIN  BZERO(I):  DIN  B0NE(1):  DIN  BTU0(1):  DIN  BTHREE(I):  DIN 
DELTRPO(I):  DIN  UELPO(I)  DIN  PRESPO(I):  OIN  T0L0(3):  DIN 
HRSSP(I):  DIN  HRSST ( 1 ) :  DIN  NELNSX(I):  DIN  TRNKSPR(I):  DIN 
TRNICDRNP(I):  DIN  PUNPSPR(I):  DIN  PUNPDRNP(I):  DIN  TQL4(1) 
360  OIN  PDENOQ(I):  DIN  PUENOO(I):  DIN  T0L4(1):  DIN  T0L5(1):  DIN 
T0L6O):  OIN  T0L7(1):  OIN  NUNOTO(I):  DIN  PHOTO(I):  DIN 
T0L3C1):  DIN  0LDCP2O):  DIN  NRXCPRES(I):  DIN  NRXPPRES(I): 
DIN  HINCPRES(I):  DIN  NIHPPRES(I) 

362  REN  Read  data  on  run. 

365  RERD  TRUEX,  FRLSEX,  THETRG,  RGRS,  T0L2,  TRUC,  NIHINC 
370  FOR  I00X-0  TO  1 

375  RERD  PNOTI(IDOX),  NUNOTI(IDOX),  T0L3(ID0X),  PHOTO(IDOX), 
NUHOTO(IDOX),  TQL4OD0X) 

380  NEXT  IDOX 


382  REfl  Read  steady  conditions  froi  the  results  file. 

385  OPEH  “I",  «1,  "Steady  Results" 

38?  REfl  Read  uar ious  parameters. 

390  INPUT  *1,  FUELS,  OXIDIZERS,  DELTRT,  NELHSS( FUELS), 

NELfISSCOX I 0 I ZERS) ,  CSTRRO,  CSTRRSLOPE,  HIXTUREO,  CTHRUST, 
LSTRR,  CHRUUOL,  GEE,  flflSSSTR,  SVSTEflTHRUST 
395  IF  NELflSS( FUELS)  >=MELflSS(QX  1 0 1  ZERS)  THEM 

flNELHS=NELHSS  ( FUELS )  ELSE  flNELI1SsNELf1SS(0X  I D I  ZERS)  'Find 
■ost  complex  system. 

400  Dlfl  LENGTH (1,nHELttS):  Din  RLPHR(1,nNELnS):  Din  NSECTIONSSd, 
HNELHX+1 ) :  Din  DISPLRCEO,  MHELJ1S) :  Din  SUEL0CITV(1,nNELnS) 
405  QS I ZES* I HT(TRUC/DELTRT)+3  'Calculate  the  size  of  the 
velocity  queue  (alloas  for  the  combustor  dead  tiae). 

410  Din  RCCEL(1,nNELnS):  Din  UELQUEUECI,  QSIZES):  Din 

SOUHD(l,nNELnS) :  DIO  D 100(1,  MELOS) :  DIO  RRER( 1 ,  HNELHS) : 
om  FR I CT ( 1 , nNELHS ) 

415  I D3X=0 

41?  REn  Read  Bore  parameters. 

420  FOR  I D0X=FUELX  TO  OXIDIZERS 

425  INPUT  «1,  RHO(IDOS),  ULLRGE(IDOS),  NPUnPS(IDOS), 

RRER  INJECT  (I  DOS),  TRUUI(IOOX),  TRUU2OD0X),  TRUINJECT(  IDOS), 
BZERO(IDOS),  BONE(IDOS),  BTUO(IDOS),  BTHREE(IDOS), 
DELTRPO(IDOS),  UELPO(IDOS),  PRESPO(IDOS) 

430  INPUT  »1,  TRNKSPR(IDOX),  TfiNKBROPC IDOS),  HfiSSTODOX), 

punpspR(ioos),  punPDRnp(iDos),  nrsspooox),  rccel(idos,o) 

435  FOR  ID1S»1  TO  NELnSX(IDOX) 

440  INPUT  *1,  RRER(  I  DOS,  ID1S),  DIRn(IDOX,  ID1S), 

LENGTH ( IDOS, ID1S),  ALPHA (IDOS, ID1S),  DISPLRCE(IDOS, ID1S), 
SUELOC I TV(  IDOS,  ID1X),  RCCELdOOS,  I01S),  SOUNDdDOS,  ID1S), 
ID2S,  FRICTdDOX,  1 01 S) 

445  NSECTIONSS( IDOS, ID1X)sID3X  'Count  up  the  totol  number  of 
sections  (for  all  elements,  i.e.  calculate  total  storage 
requirements). 

450  I03X-I03X+ID2X+1 
455  NEXT  I01X 

460  NSECTI0NSX(ID0X,f01XMD3X:  ID3X-ID3X+! 

465  NEXT  IDOS 

470  Oin  PRESSURE (I Q3S-1 ) :  Din  UELOC I TV( I D3S-1 ) 

475  FOR  IDOXsO  TO  I03X-1 

480  INPUT  <1,  PRESSURE (IDOS),  UELOC I  TV ( IDOS)  'Input  steady 


state  values. 

85  NEXT  I  DOS 
90  CLOSE  *1 

92  REH  Initialize  velocity  queue. 

95  FOR  IDOX-FUELX  TO  0XI0I2ERX 
00  FOR  I 01 X-0  TO  I HT ( T RUC/DELTRT } +3 
05  UELQUEUEODOX, ID1X)3 

UELOC I TV(NSECT I OHSX( I DOX,HELHSX( ID0X)+1)) 

10  HEXT  1 01 X 
>15  HEXT  I00X 

>1?  REH  Initialize  various  variables  needed  in  program. 
i20  FOR  IOOX-FUELX  TO  0XIDI2ERX 

>25  ELEHEHTX=HPUHPX( ID0X)-1 :  I D1 X=HSECT I ONSX ( I DOX , ELEHEHTX+ 1 ) - 1 
>30  GOSUB  630:  OLOUEL ( I OOX ) =UELOC I T V ( 1 0 1 X- 1 ) : 

OLOPRES ( I OOX ) =PRESSURE ( 1 0 1 X- 1 ) :  GOSUB  580  'Get  CP. 

>35  OLDCP 1  ( I  OOX )  “CP :  ELEHEHTXHELEHEHTX+l :  ID1X=ID1X+1:  GOSUB 
630:  GOSUB  605  'Get  Cfl. 

>40  oLocm  ( iDox)=cn 

>45  ELEHEHTX=HELHSX(IDOX):  1 0 1 X=HSECT I OHSX < I DOX , ELEHEHTX+ 1 ) - 1 
>50  GOSUB  630:  0LDUEL(I00X)=UEL0CITY(I01X-1): 

0L0PRES( IDOX)sPRESSURE( 101X-1 ) :  GOSUB  580  'Get  Cp. 

555  0LDCP2(I00X)=CP 
570  HEXT  I OOX 
575  RETURH 

580  REH  Find  CP  based  on  local  conditions. 

585  UELR®(UELOCITY(IOIX)-PSI*(UELOCITY( ID1X)-QLDUEL( IDQX)))/ 

( 1 +THETR* (UELOC I TY ( 1 0 1 X ) -OLOUEL ( I DOX ) ) )  'Caculate  Ur. 

590  PRESR3PRESSURE(I01X)-(UELR*THETR+PSI )*(PRESSURE( (DfX)- 
OLDPRES ( I DOX ) )  'Calculate  Pr. 

595  CP«PRESR+RHO< IDOX)*UELR*(SOUHD( I0QX, ELEHEHTX)* 

( RCCEL ( I DOX , ELEHEHT X ) +RCCEL ( I OOX ,  0 ) ) *DELT RT * 

S I H(RLPHR( I DOX, ELEHEHTX) )*SQUHO( I DOX , ELEHEHTX) /UELR- 
SQUHO( I DOX, ELEHEHTX)  *FR I CT( I DOX , ELEOEHTX ) *DELTRT *RBS ( UELR ) / 
( 2*0 1 flH( I DOX , ELEHEHTX ) ) )  'Calculate  Cp. 

500  RETURH 
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605  REtl  Get  CH. 

61 0  UELS»(UEL0C I  TV(  1 01  X)-PS  I  *(UELOC I TY( I 01 X)- 
UELOC I  TY(  1 01 X+  1 ) ) )/( 1  -THETR*(UELOC  ITY(IDIX)- 
UEL0CITY(ID1X+1)))  'Calculate  Us. 

615  PRESS*PRESSURE( ID1X)+(UELS*THETR-PSI )* 

(PRESSURE( I 01 X)-PRESSURE( 1 01 X+1 ) )  'Calculate  Ps. 

620  CN»PRESS-RHO(  I  D0X)*UELS*(S0UHD(  I DOX,  ELEI1ENTX) 

+(RCCEL( IDQX,ELEHENTX)+ACCEL( ID0X,0))*DELTAT* 

S I N(ALPHR( I OOX , ELEHEHTX ) ) *SOUHD ( I DOX, ELENENTX )/UELS- 
SOUHD( I DOX, ELEHENTX) *FR I CT( I DOX, ELEHEHTX)  *DELTRT*RBS( UELS ) / 
(2*D1AN( IDOX, ELENENTX)))  'Calculate  Ce. 

625  RETURN 

630  REtl  Get  PSI  and  THETR. 

635  THETR*OELTRT/(LENGTH(  I  DOX,  ELENENTX) 

/(NSECT I ONSX( I DOX, ELEHEHTX+1 )-1-  NSECT I ONSX( I DOX, ELENENTX) ) ) 

640  PS I “SOUND ( I DOX, ELENENTX )*THET A 

645  RETURN 

Boundary  Conditions.  This  section  of  code  deserves  special  attention 
as  it  applies  the  boundary  conditions  to  the  interfaces  between  system 
elements.  In  practice,  each  boundary  condition  must  be  evaluated  and  the 
resulting  equations  (some  will  be  differential  equations)  must  be  solved 
before  code  can  be  written.  A  group  of  ‘typical*  boundary  conditions  is  given 
here,  although  It  is  recognized  that  many  more  will  be  developed  by  future 
users  in  order  to  solve  more  complex  conditions. 

650  REN  Boundary  Condi  tons  Fuel  and  Oxidizer  Zero  (ullage 
pressures). 

655  ELENENTX- 1 :  GOSUB  630  'Get  Psi,  Theta. 

660  ID1X-NSECTI0NSX( 1001,1):  GOSUB  605  'Get  Ce. 

665  0LDPRES( I OOX) -PRESSURE ( ID1X) :  0LDUEL( I00X)-UEL0CITV( ID1X) 

670  PRESSURE ( 1 01 X)-ULLAGE( I OOX)  ' Nat ch  pressure . 

675  UEL0CITY( ID1X)-(PRESSURE{ ID1X)-CN)/(S0UND( I00X, ELENENTX)* 
RHO(IDOX))  'Calculate  velocity. 

680  RETURN 

685  REN  Boundary  Conditions  Fuel  and  Oxidizer  1  (area  change 


eith  relative  eotion). 

690  GOSUB  580  'Get  Cp. 

695  ID1X*ID1X+1 :  ELENENTX*ELEnENTX+1 :  GOSUB  630:  GOSUB  605 
'Get  Ce. 

700  0LDPRES( I DOX)*PRESSURE( I 01 X) :  OLDUELC I D0X)»UEL0C ITV(IDIX) 

705  BETA- AREA (  I BOX,  ELEflEHTX) /AREA  (  I DOX,  ELEflEHTX- 1 )  '  Area 

ratio. 

707  REN  Coepute  quadratic  coefficients. 

710  fl«1-BETflA2 

715  B-2*(S0UHD( IDOX,ELEHENTX)+BETR*SOUND( ID0X,ELENENTX-1 )♦ 
BETR*SUELOC  I  TV ( I  BOX ,  1  )*(S  I  H(RLPHR(  I DOX,  ELEflEHTX- 1 )  )- 
BETR*S I H( ALPHA ( I DOX , ELENEHTX) ) ) ) 

720  C«2*(Cf1-CP)/RH0( I DOX)-SUELOC I TV( I DOX, 1 )A2*(S I H( ALPHA ( I DOX, 
ELEHENTX-1 ) )-BETR*S I H(RLPHA( i DOX, ELENENTX) ) ) A2-2*S0UND( I DOX, 
ELEflEHTX- 1 )*SUELOCITV( IDOX, 1 )*(SIH(ALPHA( I OOX, ELEOEHTX- 1 ) )- 
BETR*S  I  N(RLPHR(  I  DOX,  ELEflEHTX) ) ) 

725  UEL0CITY(ID1X)«(-B+(BA2-4*fl*C)A.5)/(2*fl): 

PRESSURE( ID1X)«CH+RH0( I00X)*S0UHD( I DOX, ELEnEHTX)*UELOC I TY( I 0 
IX)  'Calculate  doenst reae  velocity  and  pressure. 

730  UELOCITYC I01X-1 )«BETfl*UELOCITY( I 01 X)-SUEL0C I TY( I DOX, 1 )* 

(SIH(ALPHA(  IOOX,ELEftEHTX-1 ))-BETA*SIH(ALPHA( I00X, ELENENTX))) 
• 

PRESSUREC I 01 X-1 )*CP-RHO( I D0X)*S0UN0( I DOX, ELEflEHTX- 1 )*UEL0C I T 
Y(IDIX-I)  'Calculate  upstreae  velocity  and  pressure. 

735  RETURN 

740  REfl  Boundary  Conditions  Fuel  and  Oxidizer  2,  5  (value). 

745  GOSUB  580  'Get  Cp. 

750  ID1X-ID1X+1:  ELEHENTX-ELENENTX+1 :  GOSUB  630:  GOSUB  605 
'Get  Ce. 

755  OLDUELC ID0X)«UEL0CITV( ID1X) :  QLOPRES(IOOX)-PRESSURECIOIX) 

760  IF  ELEHENTX-3  THEN  TRU-TRUUI(IDOX)  ELSE  TRU-TAUU2OD0X) 

765  BETR-ARER(  I  DOX,  ELEflENTX)/RRER(  I  DOX,  ELEflEHTX- 1 )  '  Area 

ratio. 

767  REN  Calculate  quadratic  coefficients. 

770  fl-1/(TRU)A2 

775  B-RHOC ID0X)*(S0UHD( ID0X,ELENEHTX-1 )*BETA+ 

SOUHOC I  DOX,  ELEflEHTX)) 

780  C-CH-CP 

785  UELOCITYC ID1X)-(-B+(BA2-4*R*C)A.5)/C2*A)  'Coepute 


‘  i 


86 


velocities. 

790  UaOCITV(  ID1X-1 ) *BETfl*UELOC I  TV ( 1 0 1 1 ) 

795  PRESSURE ( ID1X-1 )»CP-RH0( I D0X)*S0UND( I DOX, ELEflEHTX-1 )* 
UELOCITY(IDIX-I)  'Coepute  pressures. 

800  PRESSURE( I 01 X)»CN+RH0( I D0X)*S0UHD( I DOX, ELEHEHTX)* 
UELOCITY(IOIX) 

805  RETURN 

810  REN  Boundary  Conditions  Fuel  and  Oxidizer  3  (Puep  eith 
relative  notion  and  lunped  conpliance  (isothernal  bubble 
expansion)). 

815  60SUB  580  'Get  Cp. 

820  ID1X-ID1X+1:  ELEH£NTX*ELEHENTX*T :  GOSUB  630:  GOSUB  605 
'Get  Cn. 

825  OLDPRES(IDOX)-PBESSUBECIDIX):  0LDUEL(ID0X)-UEL0CITV(ID1X) 

830  BETR-BRERC I  DOX, ELENENT X ) /RREfl ( I DOX , ELENENT X- 1 ) 

835  OELTRP-DELTRPOC I D0X)-BZER0( I DOS)* 

(UELOC I TY(NSECT I 0HSX( i DOX , NPUHPX( ID0X)))-UELP0( I00X))A2- 
(B0NE( I DOX ) - ( PRESSURE ( NSECT I ONSXC i DOX, NPUNPX( I DOX) )-1 )- 
PRESPOC I DOX) )/(BTUO( I DOX)*  BTHREEC i DOX)* 

(PRESSURECNSECT I ONSXC I DOX, HPUNPXC iOQX))-1 )-PRESPO( I DOX) ) ) ) 
'Coepute  puep  pressure  rise. 

840  IF  PNOTOC I0QX)*NUN0T0( ID0X)>1  THEN  870  'Branch  for  non 
zero  conpliance. 

845  UELOC I TV( I D1 X)*(DELTRP+CP-CH+RHQ( I DOX)* 

SOUNDC I DOX, ELEHEHTX-1 )*SUEL0C I TV ( I DOX, NPUNPXC I DOX)  )* 

(SI  N( RLPHRC I DOX , ELEHEHTX- 1 ) ) -BETR*S i N(RLHR( I DOX, ELENENTX) ) ) ) 
/(RH0( ID0X)*(S0UND( IDOX,ELENENTX)+SOUHD( ID0X,ELENENTX-1 )* 
BETfi))  'Coepute  zero  conpliance  doanstrean  velocity. 

850  PRESSUREC ID1X)-CN+RH0( ID0X)*S0UHD( I DOX, ELENENTX)* 
UELOCITV(IDIX)  'Coepute  zero  conpliance  doanstrean 
velocity. 

855  UELOCITYC I01X-1 )-UEL0CITV( I01X)*BETR- 

SUELOC I TV( I DOX, NPUNPXC I DOX) )*(S I NCRLPHRC I DOX, ELEHENTX-1 ) )- 
B£TR*S I NCRLPHRC I DOX, ELENENTX) ) )  'Coepute  zero  conpliance 
upstrean  velocity. 

860  PRESSUREC ID1 X-1 )«CP-RH0( I D0X)*S0UHD( I DOX, ELENENTX-1 )* 
UELOCITV(IDIX-I)  'Coepute  zero  conpliance  upstrean 
pressure. 

865  GOTO  925  'Branch. 


870  CP!1I0*(0L0CP1  ( IOOX)+CP)/2:  Cf1HID»(0LQef11(ID0X)+CH)/2 
'Coapute  averaged  conditions. 

875  ENDPR£S*PRESSURE( ID1X-1 )  ‘First  guess. 

880  EHDPRESI*  EHDPRES:  GOSUB  1460  ‘Calculate  HendPres2. 

885  HEHDPRES 1 *HEHDPRES2 :  EHDPRES-EHDPRES*! . 00001 : 

EHDPRES2-EHDPRES:  GOSUB  1460  ‘Calculate  Hendpres2  for 
second  guess. 

890  fl*(HEHDPRES1 -HEHDPRES2)/(EHDPRES1 -EHDPRES2) : 

B*HEHDPRES1 -fl*EM0PRES1  ‘Coapute  I  inear  regression 
coefficients. 

895  EHDPRES*B/(1-fl):  EHDPRESI -EHDPRES2:  HEH0PRES1-HEHDPRES2: 
GOSUB  1460  ‘Linearly  interpolate  next  guess  and  try  it. 

900  IF  RBS(HEHDPRES2-EHDPRES) >  T0L4(I00X)  THEN  EH0PRES2-EHDPRES: 
GOTO  890  ‘Iterate  untill  sufficient  accuracy  is  achieued. 

915  PRESSURE(ID1X-1)*EHDPRES:  UEL0C I TV( I 01 X-1 )- I HUEL  ‘Coapute 

upstreaa  pressure  and  velocity. 

920  PRESSURE( I 01 X) -ENDPRES+DELTRP :  UEL0C I TV( I D1 X)=0UTUEL 
‘Coapute  dosnstreaa  pressure  and  velocity. 

925  0LDCni(ID0X)-Ct1:  0L0CP1 ( I 00X)-CP  ‘Save  old  values  for  next 

tiae. 

940  RETURH 

945  REH  Old  Chaaber  Uelocities  (Chaaber  dynaaics). 

950  OLDTHRUST-SVSTEHTHRUST  ‘Save  old  thrust. 

955  FOR  I00X-FUELX  TO  0XIDI2ERX 

960  BRCtCTIHE-TRUC:  GOSUB  1155  ‘Get  delayed  injection 
velocities  froa  queue. 

965  0L0UEL1  ( I  D0X)a0LDUEL :  BRCICT I HE-TRUC+DELTRT :  GOSUB  1155 
‘Get  older  delayed  velocity  froa  queue. 

970  0LDUEL1 ( ID0X)-(0L0UEL+0L0UEL1 ( IDQX))/2  ‘Coapute  averaged 
velocities. 

975  HEXT  ID0X 

980  CSTRR*CSTRR0+CSTRRSL0PE*(RH0 ( OX I 0 I ZERX)* 

RRER I HJECT(0X I 0 I ZERX)*0L0UEL 1(0X101 2ERX)/(RH0(FUELX)* 
ftREft INJECT (FUELX)*0LDUEL1 (FUELX) )-  HIXTUREO)  ‘Coapute 

aedicm  chaaber  exit  velocity. 

985  TCHRn«LSTRR*CSTRR/(THETRG*RGRS)  ‘Coapute  aedian  chaaber 
temperature. 


990  SVSTEHTHRUST «(RH0(0X I D I 2ERX)*RRER i HJECT (OX I D I ZERX)* 

0LDUEL1 (OX I D I ZERX ) +RHO ( FUELX ) *flREfl INJECT (FUEL!)* 

0LDUEL1 (FUELX) )*CSTRR*CTHRUST  1 Compute  eedian  thrust. 

995  1 01  X*HSECT  I OHSX(FUELX, HELHSX(FUELX)+1 ) : 

I D2X-HSECT I QHSX ( OX I 0 I ZERX , NELHSX ( OX I D I ZERX ) + 1 )  ‘Exit  plane 

indicias. 

1000  0 I DPRES*PRESSURE( I D1 I)  'Guess  eedian  pressure. 

1005  EHDPRES»PRESSURE( I 01 X)+DELTflT*(RGflS*TCHflH*(RHO(FUELX)* 

RRER i HJECT ( FUELS ) *OLOUEL 1 ( FUELX )+RH0( OX I 0 I ZERX)* 

RRER i HJECT ( OX i 0 I ZERX ) *OLDUEL 1 ( OX I 0 I ZERX ) ) /CHRHUOL- 
(RGfiS*TCHflf1*H  I  OPRES/(LSTRR*CSTRR  )  )  )  *  I  ntegrate 

endpressure. 

1010  IF  RBS ( H I DPRES- ( EHDPRES+PRESSURE ( I 01 X) )/2) >T0L2  THEH 
H I DPRES* ( EHDPRES+PRESSURE ( 101X))/2:  GOTO  1 005  ' Iterate 

untill  sufficient  accuracy  is  obtained. 

1 01 5  OLDCPRES ( FUELX ) -PRESSURE ( I 01 X) : 

OLDCPRESCOX I DI2ERX) -PRESSURE ( ID2X)  'Saue  old  ualues. 

1020  PRESSUREC I 01 X)*EHDPRES :  PRESSURE ( I D2X)-EHDPRES  'Saue  nee 

chaeber  pressure. 

1025  RETURH 

1030  REH  Chaeber  B.C's  (injector  eith  lueped  coepliance). 

1035  FOR  IDOX-FUELX  TO  0XI01ZERX 

1 040  ELEHEHTX-HELHSX(  I  BOX) 

1045  i 01 X-HSECT I 0HSX( I DOX, HELHSX( I D0X)+ 1 )-1 

1050  GOSUB  630:  GOSUB  580  ‘Get  Cp. 

1055  IF  PH0T1 ( I00X)*HUH0T1 ( I DOX) >1  THEH  1080  'Branch  for  non 
zero  coepliance. 

1057  REH  Coapute  quadratic  coefficients. 

1 060  fl-1 :  B-BH0( I D0X)*TAUI HJECT ( i DOX) A2*S0UH0( I DOX, ELEHEHTX)* 

RRER I HJECT ( I DOX)/  RREfl( I DOX, ELEHEHTX) : 

C--TRUI HJECT ( I DOX) A2*(CP-PRESSURE( I D1 X+1 ) ) 

1065  UELOC i TV( I D1 X+1 )«(-B+(BA2-4*R*C) A . 5)/(2*R) : 

UELOC i TV ( I 0 1 X ) -UELOC I TV ( I D 1 X+ 1 ) *RREfl I HJECT ( I DOX ) / 

RRER (I DOX, ELEHEHTX)  'Coepute  zero  coepliance  ueiocities. 

1070  PRESSURE* ID1X)-CP-RH0( ID0X)*S0UHD( IDOX, ELEHEHTX)* 

UELOC I TV( I D1 X)  'Coepute  zero  coepliance  upstreae  pressure. 

1075  GOTO  1125  'Branch. 

1080  CPHID-(CP+0LDCP2( ID0X))/2  'Coepute  aueage  uaiue  for  Cp. 


1085  ENDPRES1-PRESSURECID1X): 

NEHDPRES1 -FNENDPRESCFNN I DPRESCENDPRES1 ) )  'First  guess  and 
result. 

1090  ENDPRES2-EH0PRES1 *1.00001: 

NEH0PRES2*FNENDPRES (FHfll OPRES ( EHDPRES2 ) )  'Second  guess  and 

result. 

1095  fl»(NENDPRES1-NENDPRES2)/(ENDPRES1-ENDPRES2) : 

B-NENDPRES1 -fl*EHDPRES 1  'Coepute  linear  regression 
coefficients. 

1100  ENDPRES»B/(1-fl):  NEHDPRES*FHENDPRES (FHUI DPRES ( ENDPRES ) ) 
'Coepute  next  guess  and  try  it. 

1105  IF  RBS(EH0PRES-NENDPRES)>T0L3( ID0X)  THEM  ENDPRES 1 »ENDPRES2 : 
NENDPRES1 »HEHDPRES2 :  ENDPRES2-ENDPRES:  NENDPRES2-NENDPRES: 
GOTO  1095  'Iterate  untill  sufficient  accuracy  is  obtained. 
1110  PRESSUREC I 01 X)*EN0PRES  'Sage  old  value. 

1115  UEL0CITV( ID1X)>(CP-PRESSURE( ID1X))/(RH0( ID0X)* 

S0UND(ID0X,ELENENTX))  'Coepute  upstreae  velocity. 

1 1 20  UEL0C i TV( ID1X+1 )-TRU I NJECT( I DQX)*( PRESSUREC I 01 X)- 

PRESSURE(ID1X+1))A.5  'Calculate  doenstreae  velocity. 

1125  G0SUB  1180  'Load  velocity  queue  eith  nee  values. 

1130  0LDCP2C ID0X)aCP  'Save  old  value. 

1145  NEXT  I D0X 
1150  RETURN 

1155  REN  Velocity  queue  retriever. 

105?  REN  Find  bracketing  indicies. 

1160  ID2X*( INTCBRCKTIHE/DELTRT-1 )  NOD  QSIZEX)+1 
1165  ID3X*( INTCBRCICTINE/DELTRT-2)  NOD  QSIZEX)+1 
1170  0L0UEL-UELQUEUECI00X, I D2X)+(BRCKT I HE/DELTRT- I NTCBRCKT I NE/ 
0ELTRT))*  (UELQUEUEC ID0X, ID3X)-UELQUEUE( ID0X, ID2X)) 

'Retrieve  and  interpolate. 

1175  RETURN 


1180  REN  Load  velqueue. 

1185  ID2X-  ( I  NT  (SVST I  HE/DELTRT -1 )  NOD  QSIZEXW  'Coepute  index. 
1190  UELQUEUEC I00X, ID2X)*UEL0CiTV( ID1X+1 )  'Load  velocity  queue. 
1195  RETURN 


1200  REN  Structural  response  (5  degree  of  freedoe 
spring-eass-dashpot  eodel). 


1205  NINCSX-DELTRT/MNINC:  HINCSX—  IHT(-NIHCSX):  IF  HINCSK1  THEN 
H I HCSX=»  1  'Coapute  nuaber  of  ainor  increaents  corresponding 

to  the  aajor  tine  step  size  (for  greater  accuracy). 

1210  FOR  ID1X-0  TO  HIHCSX 

1215  REACT-Q  'Sub  all  reactions  (relative  to  the  uehicle 
structure). 

1220  FOR  IDOX-FUELX  TO  0XIDI2ERX 

1 225  RERCT -REACT +TRHKSPR ( I D0X)*D I SPLRCE( ID0X,1 )*TRHKDflnP( I DOX)* 
SUELOC I  TV  ( I  DOX,  1  )*PUflPSPR(  I  D0X)*D  I SPLRCE  (1 00X,  NPUriPX(  I  DOX))* 
PUI1PDRnP(  IDOX)*SUELOCITV(  ID0XIHPUI1PX(  I  DOX))  'Sub  tank  and 
puap  reaction  forces. 

1230  NEXT  I00X 

1 235  HCCEL(0, 0 )« (RERCT +OLDTHRUST + ( SVSTE11THRUST -OLDTHRUST ) * 

( I 01 X/H I HCSX)+FHEXC I TE(SVST I HE-DELTRT +QELTRT* ( ID1X/HI HCSX) ) ) 
/RRSSSTR  'Coapute  structure  acceleration. 

1240  FOR  IDQX-FUELX  TO  OXIDIZERX 

1 245  RCCEL( I DOX, 1 )--(TRHKSPR( I D0X)*0 I SPLRCE ( I DOX, 1 )♦ 

TRHKDflnP(  I  DOX) ‘SUELOC  ITV(IOOXJ)  )/t«)SST(  I  DOX) -ACCEL  (0,0) 
'Coapute  tank  accelerations. 

1250  RCCEL(ID0X,HPUnPX(100X))«-(PUnPSPR(ID0X)« 

0 1  SPLRCE ( I  DOX,  HPUt1PX(  I  DOX)  )+PUnPDflfIP(  I  DOX)* 
SUEL0CITV(ID0X,HPUnPX(ID0X)))/nflSSP(ID0X)-flCCEL(0,0) 

'Coapute  puap  accelerations. 

1255  IF  ID1X-HIHCSX  THEN  1280  'If  last  eleaent  then  branch. 

1260  DISPLRCEODOX, 1 )«OISPLflCE( I00X, 1 )+SUEL0CITV(ID0X, 1 )*DELTRT/ 

H I HCSX*RCCEL( I DOX, 1 )*(DELTRT/N I NCSX) A2/2  'Calculate 

displaceaent. 

1265  SUELOC I TV( I DOX, 1 ) -SUELOC I  TV ( I DOX, 1 )+RCCEL( IDOX, 1 )*DELTRT/ 
HIHCSX  'Calculate  velocity. 

1 270  0 1  SPLRCE(  1 00X,  HPUflPX(  1 00X)  )-D  I  SPLRCE ( I  DOX,  NPUf1PX(  I  DOX)  )♦ 
SUELOC I TV ( IOOX,NPUnPX( IDOX))*DELTRT/NIHCSX+ 

RCCEL ( I 00X , HPUOPX ( I 00X ) ) *  (DELTRT/H I HCSX ) A2/2  'Calculate 

displaceaent. 

1275  SUELOC I TV( I 00X, NPUf1PX( I DOX) )-SUEL0C I TV( I DOX, HPUOPX (IDOX))* 
RCCEL ( IDOX, NPUHPX( I DOX) )*OELTRT/H I HCSX  ' Ca I cu I ate 

velocity. 

1280  NEXT  I00X 

1285  HEXT  I 01 X 

1290  RCCELO ,0)-ACCEL(Q,0)  'Avoid  a  possible  confusion. 

1295  RETURH 


1300  REtl  Output  results. 

1305  PRINT  *1,  SVSTINE; :  PRINT  USING  '+*.*« — 

SVST I  TIE;  RCCELCO,  0) ; 

1310  PRINT  <2,  RCCEL(Q,Q); 

1315  OPEN  "fl",  *3,  "Fuel  Tank  DIsp.MO:  PRINT  *3, 

DISPLACE (FUEL!, 1 ); :  CLOSE  *3 
1320  OPEN  "A",  *3,  "Ox.  Tank  Disp.",10:  PRINT  «3, 
DISPLACE(0XIDIZER!,1); :  CLOSE  *3 
1325  OPEN  "fl",  *3,  "Fuel  Puep  Disp.",10:  PRINT  *3, 

DISPLACE (FUEL!,  NPU11P!(FUEL!)); :  CLOSE  *3 
1330  OPEN  "fl",  *3,  "Ox.  Puep  Dlsp.",10:  PRINT  «3, 

D I  SPLRCE(OX 1 0 1 ZERI,  NPUHP!(OX 1 0 1 2ER!) ) ; :  CLOSE  *3 
1335  FOR  ID0!*1  TO  NDAT! 

1340  OPEN  "A",  *!00!*2+1 ,  DATFILE$( 100!)+". PRES",  10:  OPEN  "A", 
*ID0!*2+2,  DRTFILE$(IDO!)+".UEL",  10 
1345  PRINT  *2*I00!+1,  PRESSURE(DATP!(IDO!)); :  PRINT  USING 

"♦* . ess - _  - ; PRESSURE (OATP! ( I DO!)  ) ; 

1350  PRINT  *2*ID0!+2,  UELOCITV(DATP!(IDOI))j :  PRINT  USING 

"♦*.*«* - _  " ;  UELQC I  TY(DATP!(  I  DO!) ) ; 

1355  CLOSE  *ID0!*2+1 :  CLOSE  *ID0!*2+2 
1360  NEXT  IDO! 

1365  PRINT 
1370  RETURN 

1375  REN  End  of  routine.  Chain  to  plot  routine. 

1380  CLOSE  *1 :  CLOSE  *2 

1400  BEEP:  BEEP:  BEEP  '(lake  user. 

1405  END 

1410  REN  Pogo  Data 

1415  DATR  >1,  0,  0.0015,  5.8867E+3,  1,  0.004,  0.0001 
1420  DATA  58.333E+6,  0,  0.1,  180.237E+3,  0,  0.01 
1422  REN  "Soft"  ORTA  58.333E+6,  50E-6,  0.1,  180.237E+3, 

283.13E-6,  0.01 

1425  DATA  58.333E+6,  0,  0.1,  326.93E+3,  0,  0.01 
1427  REN  "Soft"  DATA  58.333E+6,  50E-6.  0.1,  326.93E+3,  283.13E-6, 
0.01 

1430  REN  Function  Initialization. 

1435  OEF  FNEXC I TE(AT I HE) -50000 1 *S I N(AT I f1E*FREQ*2*3 .1416) 
'Excitation. 


1 440  REM  50000 ! *S I N( AT 1 HE*50*2*3 .1416) 

1 445  DEF  FHEHDPRES01 1 DPRES) -PRESSURE ( I  D1X)+DELTRT*H 1 DPRESA2/ 

(PN0T1  ( ID0X)*HUH0T1  ( IDQX))*((CPHID-HIDPRES)/(RHO( ID0X)*SQUND 
( I DOX,  ELEHEHTX)  )*RRER(  I DOX , NELNSX(  I DOX) ) -TflU I HJECT  ( I OOX)* 

01 1 DPRES-( PRESSURE ( ID1X+1 )+OLDCPRES( I DOX) )/2) A . 5* 
RRERIHJECT(IDOX))  'Used  for  injector  b.c. 

1450  DEF  FNI1 1 DPRES ( ENDPRES ) - ( EHDPRES+PRESSURE ( 1 0 1 X ) ) /2  'Used 

for  injector  b.c. 

1455  RETURN 

1460  REN  NeePres  routine. 

1465  HIDPRES*(PRESSURE( IQ1X-1 )+ENDPRES)/2  'Rueroge  pressure. 

1470  INUEL*(CPHID-HIDPRE$)/(RHO( ID0X)*S0UHD( I DOX, ELENENTX- 1 )) 

' Infloe  rate. 

1 475  0UTUEL*(H I OPRES+DELTflP-CNH I D)/(RH0( I DOX)* 

S0UND( I DOX, ELENENTX) )  'Out floe  rate. 

1 480  STORATE-ARER( I DOX, ELENENTX- 1 )*( I NUEL-BET R*0UTUEL+ 

SUELOC I TV( I DOX, NPUNPK I DOX) )*(S I N(RLPHR( I DOX, ELENENTX- 1 ))- 
BETR*SIN(ALPHA( IDOX,  ELENENTX))))  'Storage  rate  corrected 
for  relatiue  velocity. 

1 485  NENDPRES2-PRESSURE ( ID1X-1 )+DELTRT*H I DPRESA2/(PN0T0( IDOX)* 
NUN0T0( I DOX) )*ST0RRTE 

1490  RETURN 


Steady  State  Routine.  The  dynamic  response  routine  requires  that 
initial  operating  conditions  and  other  parameters  be  available  in  the  text  file 
’Steady  Results’.  The  following  program  reads  in  this  data,  computes  initial 
conditions,  and  writes  this  file. 

10  REM  Steady  State  Initialization  Prograa. 

20  GOSUB  830  'Read  data  and  initialize  variables. 

30  GOSUB  690  'Input  user  data,  etc. 

40  GOSUB  400  'Initialize  functions. 

60  GOSUB  1170  'Get  Nee  NOOT's. 

70  FOR  IDOX-FUELX  TO  OXIDIZERX:  INTND0T(I00X)-ND0T(ID0X): 

I HTHI1D0T ( I DOX ) -HI1D0T ( I DOX ) :  N00T(I00X)-1 .01*ND0T(I00X):  NEXT 
IDOX  'Save  old  values. 

00  REtl  Begin  loop. 

90  GOSUB  1170  'Get  Nee  NOOT's. 

100  IF  RBS ( NNOOT ( FUELX ) -NDOT ( FUELX ) ) < T0L0( FUELX )  RND 

RBS(NND0T(0X I D I ZERX)-N00T (OX 1 0 1 2ERX) )<T0L0(0X 1 0 1 2ERX)  THEN 
180  'Branch  is  sufficient  accuracy  has  been  obtained. 

105  REN  Coepute  next  guess. 

110  FOR  IDOX-FUELX  TO  OXIOIZERX 

115  REN  Deteraine  Linear  Regression  Coefficients. 

1 20  R-( I NTNNDOT ( I D0X)-NND0T( I DOX) )/( INTHD0T( I DOX) -ROOT ( I DOX) ) 

1 30  B- I NTNNDOT ( I D0X)-fl*( I NTNDOT ( I DOX) ) 

140  I NTNNDOT ( I 00X)-NN00T ( I 00X) :  I NTNDOT ( I DOX) -NDOT ( I 00X)  'Save 

old  guess. 

150  H00T(ID0X)-B/(1-fl)  'Linearly  interpolate  nee  guess. 

160  NEXT  IDOX 

170  GOTO  80  'Repeat. 

180  GOSUB  590  'Fill  data  registers. 

190  OPEN  "Steady  Results"  FOR  OUTPUT  RS  «1  'Open  output  file. 

195  REN  Urite  parameters  to  output  file. 

200  PRINT  *1,  FUELX;  OXIDIZERX;  DELTRT;  HELNSX(FUELX); 

NELHSX(OXIDIZERX);  CSTflRO;  CSTRRSLOPE;  HIXTUREO;  CTHRUST; 

LSTRR;  CHRHUOL;  GEE;  STRNRSS;  SVSTEHTHRUST ; 

210  FOR  IDOX-FUELX  TO  OXIDIZERX 

220  PRINT  *1,  RHO(IDOX);  ULLRGE(IDOX);  NPUHPX(IDOX); 

AREA INJECT (I DOX);  TRUUI(tDOX);  TRUU2(I00X);  TRU INJECT (IDOX); 
230  PRINT  *1,  BZER0( I DOX) ;  B0NE( I DOX) ;  BTU0< I DOX) ;  BTHREE( IDOX); 
DELTRPO(IOOX);  UELPO(IDOX);  PRESPO(IDOX);  TRNKSPR(IDOX); 
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TANKDRNPdDOX);  TflHKflRSSC  IOOZ);  PURPSPRC I  DOS) ; 
PUNPDRNPdDOX);  PUftPnfiSS(IDOX);  ACCEL( ID0X,0) 

240  FOR  1 01  Xs  1  TO  NELtlSX(IDOX) 

250  IF  ID1X»1  THEN 

D I SPLRCEC I DOX,  1 01  X)=-flCCEL(  I DOX, 0)*TflNKNRSS(  I DOX)/TRMICSPR(  I D 
OX)  ‘Compute  tank  displacements. 

260  IF  ID1X34  THEN 

0 I SPLRCE( I DOX, I D1 X)»-flCCEL( I DOX, 0)*PUNPNRSS( I D0X)/PUI1PSPR( I D 
OX)  'Coapute  pusp  displacements. 

270  PRINT  *1,  ARE  Ad  DOX,  1 01 X) ;  DIRHdDOX,  ID1X); 

LENGTH (I DOX, 1 01 X) ;  RLPHfldDOX, ID1X);  DISPLRCE( IDOX, ID1X); 
SUEL0CITV(ID0X,  ID1X);  RCCELdDOX,  ID1X);  SOUNDdDQX,  ID1X); 
NSECT I ONSX( I DOX, I D1 X) ;  FRICK IDOX, ID1X); 

280  NEXT  ID1X 

290  NEXT  IDOX 

300  FOR  IDOX>FUELX  TO  0XI0I2ERX 
310  FOR  101 X*1  TO  NELIISX(IOOX) 

320  FOR  ID2Xs0  TO  NSECT I ONSX( IDOX, ID1X) 

330  PRINT  *1,  PRESSURE( IDOX, I01X, ID2X); 

UELOCITV(IDOX, ID1X, ID2X);  'Unite  initial  conditions. 

340  HEXT  I02X 

350  NEXT  ID1X 

360  PRINT  •!,  CPRES;  NDOTdDOX)/(RRER INJECT (I DOX )*RH0(I DOX)); 

'Unite  initial  chaeben  pnessune  and  injecton  velocities. 

370  NEXT  IDOX 
380  CLOSE  *1 

390  CHAIN  "Pogo”  'Run  the  dynamic  nesponse  noutine. 

400  REN  Initial ize  al I  functions. 

410  OEF  FNDELP( IDOX, ID1X)*RH0( !D0X)*ACCEL(0,0)* 

LENGTH ( I DOX, I D1 X)*S I N(RLPHR( I DOX, I D1 X) )-RHO( IDOX)* 
LENGTHdDOX,  ID1X)*FRICT(ID0X,  ID1X)/(2*DIRN(  IDOX,  ID1X))* 
(NDOT( mOX)/(RHO(IDOX)*RRER(IOOX, I01X)))A2  'Coepute 
pnessune  nise  accnoss  eieeent. 

420  OEF  FNUEL( I DOX, I 01 X)«N00T ( I DOX )/(RH0( I DQX)*RRER( I DOX, I D 1 X) ) 
‘Local  velocity. 

430  RETURN 

440  REN  Calculate  Injecton  Pnessunes. 

450  FOR  I00X-FUELX  TO  OXIDIZERX 


455  REfl  Uork  through  eleeents  to  deters ine  injector  pressure  for 
a  particular  floe  condition. 

460  PRESSUREC 1001, 1 ,0)=ULLRGE( IDOX)  *UI lage  pressure. 

470  PRESSUREC I DO*, ! , HSECT 1 0HSX( 1 00*,  I ) )= PRESSUREC I DO*, 1,0)+ 
FHDELPC ID0X, 1 )  'Suep  pressures. 

480  PRESSUREC I  DO*, 2, 0)=PRESSURE( I DO*, 1 , HSECT I OHSXC 1 00X, 1 ) )+ 

RHOC 1D0X)*CFNUEL( IDOX, 1 )*2-FHUELC ID0*,2)A2)/2  * Area  change 

boundary  condition. 

490  PRESSUREC I DO*, 2, HSECT I OHSXC I  DO*, 2) )=PRESSUREC i DO*, 2, 0)+ 
FHDELPC I  DO*, 2)  'Pressure  at  feed  line  ualue. 

500  PRESSUREC I DO*, 3, 0)-PRESSURE( I DO*, 2, HSECT I OHSXC I DOX, 2) )- 
CFHUELC ID0*,3)/TRUU1 ( ID0X))A2  'Ualue  b.c. 

510  PRESSURE C I DO* , 3 , HSECT I OHS* C I DO* , 3 ) ) -PRESSURE  C I DOX , 3 , 0 ) + 
FHDELPC IDOX, 3)  'Puap  inlet  pressure. 

520  OELTRP-DELTRPOC ID0X)-B2ER0C I00*)*CFHUELC I DOX, 4)- 
UELP0CID0X))A2- 

(BONEC I DO*) -(PRESSUREC I DO*, 3, HSECT I OHSXC I DO*, 3) )- 
PRESPOC ID0X))/(BTU0( IDO*)+BTHREE( IDO*)* 

C PRESSUREC I DO*, 3, HSECT I OHSXC I DOX, 3) )-  PRESPOC IDO*)))) 

'Puap  pressure  rise. 

530  PRESSUREC I DO*, 4, 0) -PRESSUREC I 00*, 3, HSECT I OHSXC I DOX, 3) )+ 
OELTRP  'Puap  ouput  pressure. 

540  PRESSUREC I DOX, 4, HSECT I OHSXC IDOX, 4) )=PRESSURE( I DO*, 4, 0)+ 
FHDELPC I  DO*, 4)  'Pressure  at  throttle  ualue. 

550  PRESSUREC I DOX, 5, 0) -PRESSUREC I DOX, 4, HSECT I OHSXC I DOX, 4) )- 
CFHUELC I D0*,5)/TRUU2( I DOX) )*2  'Ualue  b.c. 

560  PRESSUREC I 00*, 5, HSECT I OHSXC I 00X, 5) )»PRESSUREC I DOX, 5, Q)+ 
FHDELPC I DO*, 5)  ' I nj  ector  pressure . 

570  HEXT  IDO* 

^  RETURH 

590  REfl  Fill  data  arrays. 

600  FOR  I00X-FUELX  TO  OXIOIZERX 

610  FOR  ID1X-1  TO  HELIISX(IDOX) 

620  FOR  ID2X-0  TO  HSECT IOHSXC IDOX, ID1X) 

630  PRESSUREC IDOX, I01X, I D2X) -PRESSURE (IDOX, ID1X,0)+ID2X/ 

HSECT I OHSXC I DOX, I DIX)* 

CPRESSURECIDOX, I 01 X, HSECT I OHSXC I DO*, ID1X))- 
PRESSUREC IDOX, I 01 X, 0) )  'Linearly  interpolate  local 

pressures. 


Spring  and  Pamper  Constanta.  The  spring  constants  of  Table  4.5  are 
calculated  to  produce  the  specified  component  natural  frequency,  wn,  which 
is  62.83  rad/sec  (10  Hz)  for  both  pumps  and  both  tanks: 


C  =  C  2VW 


[A. 15] 


The  damping  constants  of  Fig.  4.5  are  calculated  to  satisfy  the  10%  of 


critical  damping  specification  ((*0.1) 


k  =  oi;  M 

n 


[A. 16] 


This  completes  the  description  of  the  bipropellant  rocket  system  required 
by  the  P060  analysis  routine  described  in  Appendix  A. 


t 


Setting  the  nominal  pump  ouput  velocities  in  [3.17]  to  their  steady  state 


values 


V  -  Hit 

p0"  4pD2 

VPo  f  =  12^896  m/sec 
VD„  n„=  15.006  m/sec 


[A. 12] 


The  nominal  pump  pressure  rise  is  equal  to  the  difference  between  the 
steady  state  pump  output  pressure  minus  the  steady  state  pump  inlet 


pressure: 


AP  =  p  -  p 

°  3,0  3,1) 

A  P0  F  =  58.153  MPa 
APOj0x=  58.006  MPa 


[A. 13] 


The  only  remaining  pump  parameters  to  be  determined  are  the  constants 
BO,  B1,  B2,  and  B3  which  are  defined  in  [3.17].  The  constant  BO  is  set  to 
50,000  Kg/m3  for  both  fuel  and  oxidizer  pumps.  This  constant  sets  the  slope 
of  the  pump  pressure  rise  roll  off  and,  in  this  case,  results  in  the  very 
broad  roll  offs  shown  in  Fig.  4.8.  For  simplicity,  B!  is  taken  as  0.0  Pa  and 
B2  as  0.5.  It  was  found  that  these  values  allow  the  sharp  pump  pressure 
rise  roll  off  with  decreasing  pump  inlet  pressure,  as  shown  in  Fig.  4.9.  The 
values  of  83  are  5’10“6  and  2.5* JO-6  1/Pa,  respectively.  These  values 
place  the  pump  cavitation  points  (the  vertical  lines  in  Fig.  4.9)  100  KPa  and 
ZOO  KPa  below  the  nominal  pump  inlet  pressures  of  the  fuel  and  oxidizer 
systems. 

Line  Thicknesses.  The  line  thicknesses  displayed  in  Table  4.5  are 
adjusted  to  achieve  the  sonic  velocities  specified  in  Table  4.3.  This  is  done 
by  solving  [4.2]  for  the  line  wall  thickness,  (: 


Y(K/pa2 


[A.  14] 
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Propellant  Pumps.  Figure  4.6A  schematically  depicts  the  fluid  lines  of 
the  bipropellant  rocket  system  .  The  pressures  at  the  interfaces  0  (ullage 
pressure)  and  5  (combustor)  are  known.  The  upstream  pressures  at 
interfaces  3  (P3U^  and  P3UOx  )  are  the  pump  inlet  pressures  while  the 

downstream  pressures  at  interfaces  3  (  P3j)j  and  P3>dox^  are  the  PumP 

output  pressures.  These  pressures  must  be  found  in  order  to  define  the 
pump's  characteristics  in  [3.17].  Noting  that  the  feed  and  throttle  valves  are 
open  and  can  be  ignored,  the  steady  state  pressure  loss  equations22  (which 
are  just  the  steady  state  case  of  the  Navier-Stokes  equations  used  in  the 
method  of  characteristics  solution  of  Section  II)  can  be  used  to  obtain  the 
pump  inlet  and  output  pressures: 


For  a  single  pipe: 


P2=  Pt+  pgvsina  -  i/2pV2 


t2  Lr 
Dt 


♦  pgvsinocu-  \n pV2  jjLifL) 

*pgv  since  L2-i/2pV2{jL2fL2 

Pw  =  mZ4KPa  p5  mss  KPa 

Similarly, 

pw=  5SJ5I  MPa  P5D0x=  5fL2!l£  nPa 


[A.  10] 


Setting  the  nominal  pump  inlet  pressures  in  [3. 17]  to  their  steady  state 

V3 1 U6S 

Prt<r=  180.24  KPa  326.93  KPa  (A. 11] 
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For  steady  state  conditions,  the  rate  of  change  of  the  combustion 
chamber  pressure  given  by  [3.18]  must  equal  zero.  Rearanging  [3.18] 

Voic  =  — =  0  )00  m3  (A. SI 

pc 

Propellant  Inlectors.  Using  the  definition  of  injector  pressure  drop 
and  the  specified  chamber  pressure 

0.4  Pr 

= - -  =  23.33  MPa  [A.  7] 

inj  0.6  - 


The  injector  diameters  can  be  calculated  using  the  following  continuity 
relation,  the  injection  velocity,  and  the  fuel  and  oxidizer  mass  flow  rates. 
Note,  the  diameters  are  rounded  to  the  nearest  millimeter. 


4  nr* 

TTV  pV^ 


DinjjF=  45  mm 
Dinj,ox=  21  mm 


[A. 8] 


Back  substituting  these  rounded  diameters  to  obtain  the  actual  injection 
velocities  (thus  the  fuel  and  oxidizer  will  not  be  injected  at  exactly  ISO  m/sec 
but  the  mass  flow  rates  will  be  accurate)  the  injector  flow  constants  can  be 
calculated: 


-3  3  1^2 

V„j,F  =  29.66 'IQ  *  (rrf/Kg) 

. .  _3  t  1  /2 

v*u,0x=  30.79-10  (m/Kg) 


[A.  9] 


Mass  Flow  Rates.  The  total  mass  flow  rate,  M  is  given  by^u*'5:>/ 


Total  =  -4-  =  1 13.36  Kg/sec  [A.  1] 

Co  CT 

Combining  the  definition  of  mixure  ratio  and  the  fact  that  the  total  mass 
flow  rate  is  the  sum  of  the  fuel  and  oxidizer  mass  flow  rates 

rtf  =  (T^)  =  mjlKg/SBC 

[A.  2] 

ri0x=  lip  HR  =  97.0?  Kg/sec 


Since  the  fuel  and  oxidizer  tank  masses  were  assumed  equal  to  the 
masses  of  the  propellants  stored: 

ri  =  ii  xb  [a. 3] 

Mt  F  =  Kg 
Mt,ox  =  25,300  Kg 

All  component  masses  are  now  known  allowing  the  calculation  of  the 
steady  state  vehicle  acceleration,  gv: 

9v  =  Ms*MT,F*l1T,0,*Mp,F*«p<0x 
=  19.51  m/sec2 

Chamber  Properties.  From  the  definition  of  the  gas  residence  time, 


8a,  for  steady  state  conditions 


injector  is  150  m/sec.  This  high  velocity  improves  the  efficiency  of  the 
propellant  mixing  process12.  These  values  were  chosen  for  simplicity  and  in 
actual  practice,  the  injector  flow  coefficients  would  be  known.  This  would 
automatical  specify  injection  velocities  (for  a  given  injector  pressure 
drop). 

Relative  Component  Motion.  The  tanks  and  pumps  are  tied  to  the 
structure  by  means  of  springs  (lumped  stiffneses)  and  dashpots  (lumped 
damping  elements)  as  depicted  in  Fig.  4.7.  For  simplicity,  all  spring  and 
damping  constants  will  be  selected  to  produce  identical  natural  frequencies 
for  the  fuel  and  oxidizer  pump  and  tank  motions.  The  natural  frequency  is 
taken  as  10  Hz  while  the  damping  will  be  10X  of  the  critical  value.  The 
natural  frequency  was  selected  to  produce  a  resonance  within  the  limited 
range  of  excitation  frequencies  that  are  examined  in  Section  V.  The  structure 
includes  all  parts  of  the  rocket  system  other  than  the  pumps  and  tanks,  i.e. 
fluid  lines,  supporting  structure  and  payload.  The  structure  mass  is  3000  Kg 
while  both  pump  masses  are  taken  as  100  Kg.  The  tank  masses  are  assumed 
equal  to  the  mass  of  the  propellants  contained,  neglecting  the  tanks' 
structural  masses.  For  this  system,  the  tanks  contain  enough  propellant  for 
200  sec  of  burn  time,  therefore,  the  tank  structural  masses  constitute  an 
insignificant  part  of  the  total  tank  masses.  Note  that  the  structural  mass 
fraction  (including  payload  as  part  of  the  structure)  is  approximately  12% 
which  is  typical  for  an  OTV  mission1 1 . 


Calculated  Yalugs 

The  bipropellant  POGO  analysis  routine  detailed  in  Appendix  A  requires 
still  more  parameter  values  (see  Tables  4.4  and  4.5).  in  this  section, 
relationships  taken  from  the  literature  will  be  used  to  calculate  these  values 
from  those  initial  values  given  to  define  the  system. 


propellants  to  obtain  high  specific  impulses  (high  exhaust  velocities).  The 
fuel  density  is  taken  as  71  Kg/m3  and  the  bulk  modulus  as  102.6  MPa  while  the 
corresponding  oxidizer  values  are  1 140  Kg/m3  and  1.642  GPa  *1 .  The  steady 
state  mixture  ratio  (oxidizer  mass  flow  rate  devlded  by  the  fuel  mass  flow 
rate)  is  6:1  (8:1  being  Stoichiometric) *  1  »21.  As  described  by  [3.19],  the 
effective  exhaust  velocity  is  modeled  as  a  linear  function  of  mixture  ratio 
whose  slope,  C*s|ope  is  taken  as  73  m/sec  while  the  gas  constant,  R  of  the 
exhaust  gas  is  5.89  K J/Kg •  K  for  purposes  of  this  analysis21.  Both  the  fuel 
and  oxidizer  ullage  pressures,  P^,  (the  gas  pressure  at  the  top  of  both 

tanks)  are  150  KPa  as  in  the  Dorsch  example1 . 

Fluid  Lines  and  Tanka.  The  lengths,  diameters,  sound  speeds, 
orientations  and  friction  factors  of  the  system  lines  and  tanks  given  in  Table 
4.3  represent  a  typical  system.  Note  that  the  feed  lines  are  longer  and 
larger  in  diameter  than  the  high  pressure  discharge  lines.  Also,  the  feed 
lines  are  oriented  parallel  to  the  rockets  longitudinal  axis  which  will  make 
the  pump  inlet  pressures  sensitive  to  vehicle  acceleration  variations,  since 
the  hydrostatic  pressure  at  the  bottom  of  the  feed  lines  (pump  inlet 
pressure)  is  proportional  to  the  vehicle  acceleration.  The  sound  speeds  in 
the  fuel  lines  are  chosen  to  be  similar  to  those  of  the  monopropellant  system 
modeled  by  Dorsch,  et  al. 1 ,  while  the  oxidizer  sonic  velocities  are  taken  as 
approximately  one-half  as  large  as  the  fuel  sonic  velocities  simply  to  provide 
a  contrast  between  the  fuel  and  oxidizer  systems.  The  sonic  velocities  are  a 
function  of  the  Young’s  modulus  and  thickness  of  the  line  material.  The  line 
material  Young's  modulus  is  10  GPa  (a  representative  value)12.  The  tank 
lengths  correspond  to  the  volumes  required  to  supply  200  seconds  of  steady 
state  propellant  flow  (see  calculated  parameters).  The  friction  factors  are 
taken  from  the  Dorsch  system1 . 

Y«lv»  and  Injectors.  For  this  analysis,  both  the  feed  and  throttle 
valves  (Fig.  4.6)  are  assumed  to  be  open.  This  is  accomplished  by  using  a 
relatively  large  number  ( 100  (m3/ Kg) 1/2  in  this  analysis)  for  the  valve  flow 
constants  of  the  P060  routine.  The  propellant  injector  pressure  drops  are 
taken  as  40*  of  their  respective  upstream  (end  of  discharge  ltne)  pressures 
(see  Dorsch  monopropellant  system1).  The  injection  velocity  for  each 


Appendix  B:  Bipropellant  Liquid  Rocket  System  Parameters 


In  this  appendix,  parameter  values  describing  the  bipropellant  liquid 
rocket  system  depicted  in  Figs.  4.6  and  4.7  will  be  calculated  from  those 
values  given  in  Table  4.3.  The  values  defining  the  system  will  first  be 
justified  and  then  liquid  rocket  relationships  will  be  given  for  the  calculation 
of  other  parameters  required  by  the  bipropellant  P060  analysis  routine 
documented  in  Appendix  A. 


Initial  Values 

Before  any  numerical  analysis  can  be  attempted,  values  must  be  assigned 
to  parameters  describing  the  system.  The  bipropellant  liquid  rocket 
described  in  this  thesis  is  intended  to  represent  a  'typical*  system.  As  such, 
its  descriptive  parameters  do  not  necessarily  correspond  to  any  actual 
rocket  and  this  system  is  examined  solely  to  demonstrate  the  application  of 
the  bipropellant  POGO  analysis  routine  developed  in  this  investigation. 

Rocket  Engine  Properties.  The  steady  state  thrust  of  the  system  is 
taken  as  500,000  N  with  a  steady  state  effective  exhaust  velocity,  C0*  of 
4,415  m/sec  (corresponding  to  a  specific  impulse  of  450  sec)  and  thrust 
coefficient,  Cj  of  0.999.  The  steady  state  chamber  pressure  is  35  MPa  (5000 

psia)  while  the  chamber  temperature  is  3500  K.  These  values  are  typical  of 
the  high  performance  OTVs  and  OMVs  being  researched  today11*16.  The 

values  for  the  chamber  dead  time,  X.  and  gas  residence  time,  0_  are  taken 

m  y 

as  4  msec  and  1.5  msec,  respectively.  These  values  are  very  similar  to 
those  for  the  monopropellant  system  analyzed  by  Dorsch,  et  al. 1 . 

Propel  1  ants.  The  system  employs  liquid  hydrogen  fuel  and  liquid 
oxygen  as  the  oxidizer  since  most  OTV  and  OMV  concepts  employ  these 


Steady  State  Data  Statements.  The  following  values  must  appear  in 
the  data  statements  in  the  order  shown:  FUELS,  OXIDIZERS,  CSTARO, 
CSTARSLOPE,  MIXTUREO,  CTHRUST,  LSTAR,  CHAM  VOL,  STRMASS. 

The  following  set  of  data  must  next  appear  twice:  first  for  the  fuel  system 
and  second  for  the  oxidizer  system:  NELMSS,  BULK, RHO, ULLAGE,  NPUMPS, 
AREA  INJECT,  TAUV1,  TAUV2,  TAUINJECT,  BZERO ,  BONE,  BTWO, 
BTHREE,  DELTAPO,  VELPO,  PRESPO,  TOLI,  TANKMASS,  TANKSPR, 
TANK  DAMP,  PUMPMAS5,  PUMPSPR,  PUMPDAMP. 

The  following  data  must  appear  once  for  each  element  in  the  fuel  and 
oxidizer  systems:  YOUNG,  THICK,  AREA,  DIArt,  LEN6TH,  ALPHA,  FRICT. 

Dynamic  Response  Data  Statements.  The  following  data  must  appear 
first:  TRUES,  FALSES,  THETA6,  R6AS,  TOLZ,  TAUC,  MININC. 

The  following  data  must  appear  once  for  the  fuel  system  and  once  for  the 
oxidizer  system:  PNOTI,  NUNOTI,  TOL3,  PNOTO,  NUNOTO,  TOLI. 


Array  Variables. 


ACCEL:  component  acceleration 
AREA:  component  area 
BONE:  pump  coefficient 
BTWO:  pump  coefficient 
BZERO:  pump  coefficient 
DELTAPO:  nom.  pump  pres,  rise 
DISPLACE:  component  displacement 
FALSE!:  logical  value  («0) 

FUELS:  fuel  index  (=0) 

INTNMDOT:  intermediate  mass  rate 
MOOT:  current  mass  rate 
NMDOT:  final  mass  rate 
NSECTIONSS:  section  index 
NUNOT1:  injector  bubble  nom.  volume 
OLDCPi:  old  Cp  value 
OLDCPRES:  old  chamber  pressure 
PNOTO:  pump  bubble  nom.  pres. 
PRESPO:  nom.  pump  inlet  pres. 
PUMPDAMP:  pump  damping  constant 
PUHPSPR:  pump  spring  constant 
SOUND:  sonic  velocity 
TANKDAMP:  tank  damping  constant 
TANKSPR:  tank  spring  constant 
TAUVi:  feed  line  valve  flow  const. 
THETAS:  gas  residence  time 
TOLf:  pump  match  pres,  tolerance 
TOLi:  pump  inlet  pres,  tolerance 
ULLAGE:  ullage  pressure 
VELPO:  nom.  pump  velocity 
YOUNG:  component  Young’s  modulus 


ALPHA:  component  orient,  angle 
AREAINJECT:  injector  area 
BTHREE:  pump  coefficient 
BULK:  fluid  bulk  modulus 
DELMSS:  element  index 
DIAM:  component  diameter 
ELEMENTS:  element  index 
FRICT:  component  friction  factor 
INTMDOT:  intermediate  mass  rate 
LENGTH:  component  length 
NELMS:  ft  of  elements  in  system 
NPUMPS:  pump  element  ID  * 
NUNOTO:  pump  bubble  nom.  volume 
OLDCMI:  old  Cp  value 
OLDCPZ:  old  Cp  value 
OXIDIZERS:  oxidizer  index  (-1 ) 
PNOTI:  injector  bubble  nom.  pres. 
PRESSURE:  pressures 
PUMPMASS:  pump  mass 
RHO:  fluid  density 
SVELOCITY:  structure  velocity 
TANKMASS:  tank  mass 
TAUINJECT:  injector  flow  constant 
TAUVZ:  throttle  valve  flow  const. 
TOLO:  mass  flow  match  tolerance 
TOL3:  injector  match  pres.  tol. 
TRUES:  logical  value  (  —  1) 
VELOCITY:  velocities 
VELQUEUE:  velocity  queue 


Program  Variables 


This  section  lists  all  of  the  variables  used  in  the  steady  state  and  dynamic 
response  routines.  The  identity  and  meaning  of  each  program  variable  will 
be  given,  and  the  input  order  will  be  specified  for  those  variables  whose 
values  are  read  from  data  statements. 


Simple  Variables. 


ATIME:  time  parameter 

BETA:  area  ratio 

CM:  current  value  of  Cm 

CP:  current  value  of  Cp 

CPRES:  chamber  pressure 

CSTARO:  nominal  exhaust  velocity 

CTHRUST:  thrust  coefficient 

DELTAT:  time  step 

ENDPRES:  pressure  (intermediate) 

ENDPRES2:  pressure  (int. ) 

ID  IX:  counter 
ID3Z:  counter 
INVEL:  pump  inlet  velocity 
MIDPRESS:  pressure  (average) 
MIXTURE:  current  mixture  ratio 
MNELMZ:  largest  t  of  elements 
NENDPRESl:  pressure  (intermediate) 
NEWPRESS:  pressure  (terminal) 
OLDTHRUST:  last  thrust  value 
PRESR:  Pr 

PSI:  grid  parameter 
REACT:  total  reaction  force 
R6A3:  chamber  gas  constant 
SIDEC  designator  (fuel-0,  ox. si) 
STRMASS:  structure  mass 
TAU:  orafice  constant 
TCHAM:  chamber  temperature 
THETA:  grid  mesh  ratio 
TOL2:  combustor  pres,  tolerance 
VELR:  Vr 


BACKTIME:  Velocity  queue  param. 
CHAMVOL:  chamber  volume 
CMMID:  Cm  (intermediate) 

CPMID:  Cp  (Intermediate) 

CSTAR:  current  exhaust  velocity 
CSTARSLOPE:  local  slope 
DELTAP:  pump  pressure  rise 
DELTATOUT:  output  time  step 
ENDPRES  I :  pressure  (int. ) 

IDOS:  counter 
ID2S:  counter 

INTPRES:  pressure  (intermediate) 
LSTAR:  chamber  char,  length 
MININC:  min.  structural  time  step 
MIXTUREZERO:  nom.  mixture  ratio 
MNSECTX:  largest  »  of  sections 
NENDPRESZ:  pressure  (Int.) 
NINCSX:  *  of  minor  timesteps 
OUTVEL:  pump  outlet  velocity 
PRESS:  Ps 

QSIZEX:  size  of  velocity  queue 
RES$:  a  user  response  string 
SECS:  section 

STORATE:  pump  inlet  storage  rate 

5YSTEMTHRUST:  thrust 

TAUC:  chamber  dead  time 

TEND:  end  time 

THICK:  pipe  wall  thlcknes 

TSART:  begin  time 

VELS:  Vs 


4 


1080  ORTA  1E+10,  644.37E-6,  70.686E-3,  0.3,  8,  1.5708,  0.03 

1090  DATA  1E+10,  644.37E-6,  70.686E-3,  0.3,  8,  1.5708,  0.03 

•  1100  DATA  1E+10,  1.9718E-3,  17.671E-3,  0.15,  1,  0,  0.002 
1110  DATA  1E+10,  1.9718E-3,  17.671E-3,  0.15,  1,  0,  0.002 
1120  DATA  1E+10,  1.9245E-3,  7.0686,  3,  2.4093,  1.5708,  0.03 
1130  DATA  1E+10,  305.05E-6,  22.698E-3,  0.17,  3,  1.5708,  0.03 

•  1140  DATA  1E+10,  305.05E-6,  22.698E-3,  0.17,  3,  1.5708,  0.03 

1150  DATA  1E+10,  506.58E-6,  5.6745E-3,  0.085,  1,  0,  0.002 

1160  DATA  1E+10,  506.58E-6,  5.6745E-3,  0.085,  1,  0,  0.002 

1170  RED  Thrust  deters i nation. 

•  1180  n I XTURE-nDOT (OX  1 0 1 ZERX) /HOOT ( FUELX )  'Cospute  sixture 

ratio. 

1185 

SVSTEflTHRUST  * ( HOOT ( FUEL! ) +HDQT ( OX i 0 1 ZERX) )*(CSTARO+CSTARSLOP 
4  E*  (ROOT (OX 1 0 1 ZERX) /HOOT ( FUELX )-fl I XTUREO) ) *CTHRUST 

' Deters ine  thrust. 

1190  ACCEL ( FUELX , 0 ) -SYSTERTHRUST/ ( STRHASS+TAHKHASS ( FUELX ) ♦ 
TAHKRASSCOX 1 0 1 ZERX ) +PUflPRflSS ( FUELX ) +PU  WR1ASS ( OX 1 0 1 ZERX) ) 

‘Get  systes  acceleration. 

•  1200  ACCEL ( OX I D 1 ZERX , 0 ) *RCCEL ( FUELX , 0 ) 

1210  CSTRR-CSTRRO+CSTRRSLOPE* ( R I XTURE-H I XTUREO )  ‘Operat ing 

exhaust  velocity. 

1 220  CPRES»LSTAR*CSTAR/CHAHUOL*(HOOT (FULEX)+RDQT (OX  1 0 1  ZERX) ) 

•  'Operating  chasber  pressure. 

1230  GOSUB  440  'Get  Injector  Pressures. 

1240  FOR  IDOX-FUELX  TO  OXIOIZERX 
1250  NRDOT ( I DOX)*TAU I HJECT ( I DOX)* 

(PRESSURE ( I DOX, HELHSX( I DOX) ,  HSECT I OHSX( I DOX, HELRSX( I DOX) ) )- 
CPRES)* . 5*AREA I HJECT ( I DOX)*RHO( I DOX)  'Cospute  boss  floe 

rate. 

1260  NEXT  I DOX 
1270  RETURH 

«  1280  EHD 
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870  FOR  IDOX»FUELX  TO  0XIDI2ERX 

880  RERO  NELHSX(IDOX),  BULK(IDOX),  RHO(IDOX),  ULLAGE(IDOX), 
NPUNPX(IQOX),  AREAINJECT( IDOX),  TAUUI(IDOX),  TflUU2(ID0X), 
TRU INJECT (I DOX),  BZERO(IDOX),  COHE(IDOX),  BTUO(IOOX), 
BTHREE(IDOX) 

890  RERO  DELTRPQ(IDOX),  UELPO(IOOX),  PRESPO(IOOX),  TOLO(IDOX), 
TRNICnflSSC  IDOX),  TANKSPR(  IDOX),  TANKDANP(  I OOX) , 

PUNPNRSS( I DOX) ,  PUNPSPR( IDOX),  PUNPDRHP( I OOX) 

900  NEXT  I OOX 

910  IF  NELNSX ( FUELX )  >NEU1SX ( OX  1 0 1 2ERX )  THEN  NNELNXsHELHSX ( FUELX ) 
ELSE  NNELHX*NELNSX ( OX I D I 2ERX ) 

920  DIN  AREA (1, NNELNX):  DIN  D I RN(  1  y MHELflX) :  DIN 
LEHGTH(  1 , flNELMX) :  DIN  ALPHA (1, NNELNX):  DIN 
0 I SPLRCEC1 , NNELNX) :  DIN  SUELOC I TY(1 , NNELNX ) :  DIN 
ACCEL (1, NNELNX):  DIN  SOUND (1 .NNELNX) :  DIN  HSECTIONSXd, 
NNELNX):  DIN  FRICTd,  NNELNX) 

930  FOR  IDOX=FUELX  TO  0XI0I2ERX 
940  FOR  I01X=1  TO  NELNSX(IDOX) 

950  READ  YOUNG,  THICK,  AREA  (IDOX,  ID1X),  DIRNdDOX,  ID1X), 
LENGTHdOOX,  101 X) ,  ALPHA  (I  DOX,  I01X),  FRICT(  IDOX,  ID1X) 

960  SOUNDdDOX, 1D1X)-(BULK( ID0X)/(RH0(ID0X)*(1+BULK( IDOX)* 
0IAHd00X,ID1X)/(Y0UNG*THICK))))\5  'Coipute  local  sonic 
velocity. 

970  NEXT  I01X 

980  n00T(ID0X)-UELPQ( IDOX)*RRER(IDOX,NPUnPX(IDOX))*RHO(IOOX) 
'First  guess  for  systee  floe. 

990  NEXT  IDOX 
1000  RETURN 

1010  REN  Data  for  steady  state  problea. 

1020  DATA  0,  1,  4415,  75,  6,  0.999,  7,  0.1,  3000 
1030  ORTA  5,  102.24E+6,  71,  150E+3,  4,  1.5904E-3,  100,  100, 
29.664E-3,  50000,  0,  0.5,  5E-6 
1040  DATA  58.153E+6,  12.896,  180.24E+3,  0.1,  3236,  12.774E+6, 
40.663E+3,  100,  394.78E+3,  1.2566E+3 
1050  ORTA  5,  1 .6416E+9,  1140,  150E+3,  4,  572.56E-6,  100,  100, 
30.787E-3,  50000,  0,  0.5,  2.2033E-6 
1060  DATA  58.006E+6,  15.006,  326.93E+3,  0.1,  19414,  76.643E+6, 
243.96E+3,  100,  394.78E+3,  1.2566E+3 
1070  ORTA  1E+10,  2.0448E-3,  7.0686,  3,  6.4478,  1.5708,  0.03 


640  UELOCITVdDOX,  ID1X,  1  D2X)«FHUEL(  1 DOX,  ID1X)  'Element 
velocities. 

650  NEXT  ID2X 
660  NEXT  ID1X 
670  NEXT  IOOX 
680  RETURN 

690  REN  Input  user  data  and  calculate  required  nuaber  of 
sect i ons . 

700  CLS:  CALL  TEXTFONT(O):  CALL  TEXTS I 2EC 1 8) :  PRINT:  PRINT 
TRB(5); "Steady  State  Initialization” 

710  CALL  TEXTSI2EC 12) :  PRINT:  PRINT 

720  PRINT  TRB(5);”Tiae  Increaent:  INPUT  ”",DELTRT 

730  HNSECTX=0 

740  FOR  IDOX-FUELX  TO  0XI012ERX 
750  FOR  ID1X*1  TO  NELNSX(IOOX) 

760  NSECTIQNSXdDOX, ID1X)=INT(LENGTH(IOOX, ID1X)/(DELTRT* 

( SOUND ( 1001,101  X)+t100T  ( 1 00X)/(RH0(  ID0X)*RRER(  I DOX, I 01 X) ) ) ) ) 
'Coapute  nuaber  of  sections  based  on  Courant  condition. 

770  IF  NSECTIQNSXdDOX,  ID1X)*0  THEN  PRINT  TRB(5);”Tiae  Increaent 
too  large.*:  END  'Quit  if  response  is  unacceptable. 

780  IF  NSECT 1 0NSX( I DOX, 1 01 X) >HHSECTX  THEN 

nNSECTX=NSECTIONSX(IDQX, 1 01 X)  'Oeteraine  systea  aith 
aaxiaua  nuaber  of  sections. 

790  NEXT  ID1X 
800  HEXT  I00X 

810  DIN  UEL0CITY(1  ,HNELHX,f1NSECTX) :  DIM 
PRESSURE  ( 1 ,  HNELflX,  flHSECTX) 

820  RETURN 

830  REH  Input  data  and  diaension  arrays. 

840  READ  FUELX,  OXIOIZERX,  CSTflRO,  CSTHRSLOPE,  MXTUREO, 

CTHRUST,  LSTRR,  CHflflUOL,  STRNRSS 

850  Din  NELnSX(l):  Din  RH0(1):  Din  ULLAGE(I):  Din  NPUnPX(l):  DIN 
RRERINJECT(I):  DM  TRUUI(I):  DM  TRUU2C1):  DM  TRUINJECT(I): 
DM  n00T(1):  DM  BZEROd):  DM  B0NE(1):  DM  BTU0(1):  DM 
BTHREE(I):  DM  INTnOOT(l):  OM  NttDOTd):  DM  I NTNODOT ( 1 ) 

860  DM  DELTAPO(I):  DM  UELPO(I):  DM  PRESPO(I):  OM  PUnPSPR(l): 
DM  PUnPORnPd):  OM  TRNXSPR(t):  DM  TRNKDflnP(  1 ) :  DM 
PUnPHRSS(  1 ) :  DM  TflNiCnnSS( I ) :  DM  BUUC( 1 ) 
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In  this  thesis,  a  generalized  computer  code  for  the 
simulation  of  POGO  instability  in  liquid  rockets  was  gener¬ 
ated  and  verified.  The  term  POGO  refers  to  the  out  of 
phase  motion  of  the  ends  of  a  liquid  rocket,  which  can 
build  to  a  dangerous  magnitude  due  to  propulsion  feed  back. 
The  elimination  of  any  POGO  instability  will  be  of  great 
importance  in  the  development  of  future  liquid  rocket  sys¬ 
tems  such  as  OTV ' s  and  OMV's.  The  model  detailed  in  this 
investigation  is  based  on  a  method  of  characteristics  solu¬ 
tion  of  the  simplified  Navier-Stokes  equations.  The  resul¬ 
ting  non-linear  differential  equations  were  solved  using 
first  order  finite  difference  methods.  This  solution  was 
applied  to  the  fluid  lines  of  several  liquid  rocket  systems. 
Boundary  conditions,  based  on  component  models  used  gener¬ 
ally  throughout  the  literature  were  developed  to  link  these 
.fluid  lines  . 

The  computer  routines  were  verified  by  comparison  with 
published  results  from  several  sources.  The  results  of  nu¬ 
merous  runs  agreed  quite  well  with  the  published  data,  even 
in  very  non-linear  systems  with  large  disturbance  amplitudes. 

As  an  aid  to  future  investigators  the  routines  developed 
in  this  thesis  were  applied  to  a  typical  bipropellant  liquid 
rocket  system.  Both  the  transient  and  the  steady  state  re¬ 
sponses  of  this  system  to  sinusoidal  thrust  variations  were 
obtained.  A  relatively  limited  parametric  study  was  perform¬ 
ed  and  indicated  the  extreme  parametric  sensitivity  of  the 
model  to  input  thrust  variation.  Two  factors  were  most  im¬ 
portant:  The  pump  operating  points  and  the  chamber  dead 

time,  the  stability  of  the  system  being  greatly  decreased 
with  a  large  dead  time  or  by  operation  near  the  cavitation 
point . 

A  user's  guide  was  compiled  which  details  the  software 
developed  here.  Its  purpose  is  to  facilitate  the  application 
of  these  routines  to  other  systems  by  future  investigators. 
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